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I.  INTRODUCTION 


A.  MOTIVATION 

Historically,  robotics  research  has  focused  on  individual  capabilities,  such  as 
traveling  between  points  and  obstacle  avoidance.  However,  a  fundamental  shift  is 
occurring:  robots  are  increasingly  being  put  to  work  in  real-world  environments.  These 
environments  tend  to  be  complex  and  cluttered,  and  the  tasks  are  complicated,  requiring 
advances  in  controls,  sensing,  perception,  and  communication.  In  particular,  dive 
operations  are  inherently  dangerous.  Physiological  effects  limit  dive  duration  and 
frequency  and  necessitate  a  large  support  crew,  increasing  operational  costs.  The  sensory- 
deprived  underwater  environment  makes  navigation,  communication,  and  documentation 
challenging.  The  Center  for  Autonomous  Vehicle  Research  (CAVR)  at  the  Naval 
Postgraduate  School  (NPS)  is  developing  a  Robotic  Diver  Assistant  System  (RDAS)  to 
provide  autonomous  support  to  diver  teams,  which  has  the  potential  to  significantly 
enhance  underwater  operations.  The  RDAS  project  is  aimed  at  providing  utility  to  the 
diver  team  (e.g.,  illumination,  improved  situational  awareness,  etc.)  without  burdening 
the  team  with  vehicle  command  and  control,  thereby  augmenting  the  diver  team  and 
allowing  more  effective,  efficient,  and  safer  operations.  This  program  seeks  to  go  beyond 
co-inhabitance  of  man  and  machine — the  aim  is  to  fundamentally  enable  the 
transformative  capability  of  robots  as  underwater  co-workers. 

The  RDAS  finds  application  in  many  naval  operations,  but  of  particular  interest  is 
the  potential  benefit  to  the  salvage,  explosive  ordinance  disposal  (EOD),  and  undersea 
rescue  operations  of  the  Department  of  the  Navy.  This  application  requires  operation  of  a 
hovering-class  unmanned  underwater  vehicle  (UUV)  in  close  proximity  to  humans  as 
well  as  other  features  (e.g.,  structures,  the  sea  bottom,  etc.),  which  in  turn  requires 
precision  control  of  the  vehicle.  The  application  will  require  testing  of  a  vehicle  in 
multiple  configurations  (i.e.,  different  payload  combinations)  as  well  as  the  development 
of  various  perception  and  control  strategies.  An  accurate  dynamic  model  is  required  to 
facilitate  precision  control.  However,  these  hydrodynamic  models  are  notoriously  time 

and  resource  consuming  to  develop  in  practice.  Developing  individual  models  for  all 
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possible  configurations  is  impractical.  An  effective  and  efficient  method  for  learning  the 
dynamic  model  online  is  necessary  to  execute  the  RDAS  program  as  well  as  other  closed- 
quarters  operations. 

The  goal  of  this  research  is  to  learn  a  parameterized  dynamic  model  for  an  UUV 
in  real  time,  and  to  automate  the  procedure  where  possible.  Each  change  in  payload 
configuration  changes  the  weight  distribution,  center  of  buoyancy,  and  drag  properties  of 
the  platfonn  (i.e.,  changes  the  dynamics),  requiring  a  new  model.  Initially,  a  propulsion 
model  is  developed  using  experimentally  obtained  thruster  data,  then  the  full  equations  of 
motion  for  a  free  body  are  developed.  Then,  several  assumptions  are  made  to  simplify 
these  equations  which  are  in  turn  used  to  create  a  vehicle  simulator.  Various  techniques 
in  System  Identification  and  probabilistic  state  estimation  (that  account  for  uncertainties 
in  the  measurement  and  modeling  process)  are  investigated  to  develop  a  parametric 
model  of  the  system  as  well  as  a  non-parametric,  neural-network-based  approach.  These 
methods  are  investigated  to  determine  which  of  the  four  are  able  to  be  used  for  this 
application  as  well  as  the  strong  and  weak  points  of  these  methods.  The  applicability  of, 
and  results  from,  these  results  are  tested  using  the  developed  simulator  and  then  either 
applied  to  the  vLBV300  or  discarded.  This  work  culminates  in  a  hydrodynamic  model  of 
the  SeaBotix  vLBV300  remotely  operated  vehicle  (ROV),  which  is  the  development 
platform  for  the  CAVR  robotic  diver  assistant. 

B.  LITERATURE  REVIEW 

An  accurate  dynamic  model  for  an  UUV  is  required  to  inform  the  parameter 
learning  process,  for  simulation  of  a  vehicle,  and  for  validation  of  experimental  results. 
Fossen  [1]  presents  a  detailed  approach  exploring  the  kinematics  and  dynamics  that  make 
up  a  complete  model  for  a  generic  six  degree  of  freedom  underwater  vehicle.  He  presents 
an  approach  using  the  Newton-Euler  formulation  based  on  Newton’s  second  law 

™vc=fc  (!) 

where  m  is  the  mass  of  a  body,  t)c  is  acceleration  and  fc  is  force  acting  on  the  body.  It 
can  easily  be  seen  from  Equation  (1)  that  if  no  force  is  acting  on  the  body  then 
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acceleration  will  be  zero  and  the  body  will  move  at  a  constant  speed.  Fossen  then 
describes  Euler’s  first  and  second  axioms  to  Newton’s  second  law, 


<11 

o 

Pc  =  mvc 

(2) 

hc  =  mc 

hc  —  Ic co 

(3) 

where  fc  and  mc  are  the  forces  and  moments  acting  on  a  body  referenced  to  its  center 

of  gravity,  Ic  is  the  inertia  tensor  about  the  body’s  center  of  gravity,  and  co  is  the  angular 
velocity  vector.  These  axioms  describe  Equation  (1)  in  terms  of  the  conservation  of 
linear,  pc,  and  angular  momentum,  hc ,  which  is  very  convenient.  Fossen  applies  these 

concepts  to  develop  his  generic,  six  degree  of  freedom  (DOF),  nonlinear,  and  dynamic 
equations  of  motion.  This  generic  model  is  then  specialized  for  the  development 
platform  to  account  for  the  unique  design,  relatively  slow  speeds  of  operation,  and  high 
level  of  maneuverability.  The  results  from  [1]  are  extensively  used  in  this  research. 

System  identification,  also  known  as  model  learning,  attempts  to  leam  model 
parameters  by  systematically  exciting  various  dynamic  modes  of  the  vehicle  and  is  an 
established  research  field.  System  identification  techniques  are  distinguished  by  when 
they  are  applied  to  gathered  data  (online  vs.  offline)  and  what  type  of  model  is  of  interest 
(parametric  vs.  non-parametric). 

Offline  approaches  capture  the  vehicle  motion  for  a  commanded  behavior  and 
then  choose  parameters  that  best  represent  all  the  data  offline  (i.e.,  regression  techniques, 
such  as  Linear  Least  Squares  Regression,  which  minimize  the  compounded  error  for  the 
whole  data  set).  Slotine  and  Li  [2]  and  Astrom  and  Wittenmark  [8]  present  various 
analytical  and  computational  methods,  including  the  gradient  estimator  (which  minimizes 
the  instantaneous  estimation  error)  and  versions  of  the  least  square  estimator  for  system 
identification  to  handle  time-varying  parameters  and  ill-conditioned  learning  problems. 

Alternatively,  model  parameters  can  be  estimated  online:  the  estimate  is 
continually  updated  as  new  data  becomes  available.  The  estimation  process  itself  is 
considered  as  a  dynamic  process,  to  which  stability  and  convergence  analyses  can  be 
applied.  It  is  often  possible  to  adapt  offline  techniques  for  online  implementation  (e.g., 
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recursive  least  squares  regression),  or  to  develop  new  techniques  entirely  (e.g.,  Gradient 
Descent  methods)  [2,  8].  Andrieu  et  al.  [12]  investigates  a  method  for  online,  point 
estimation  of  static  parameters  for  general  state-space  models.  Kugler  [13]  presents  a 
non-linear  parameter  to  output  operator  approach  that  allows  it  to  be  applied  to  both  finite 
and  infinite  dimensional  differential  equations.  Thrun  and  Bulgard  [3]  present  a 
probabilistic  approach  based  on  the  Bayes  filtering  process  for  online  state  estimation  that 
accounts  for  noise  and  uncertainty  in  the  state  and  process  measurements.  The  Bayesian 
approach  is  premised  on  the  concepts  of  controllability  and  observability.  A  method  for 
determining  both  the  controllability  and  observability  of  a  system  is  presented  by  Ogata 
[18]. 

Offline  system  identification  techniques  have  been  applied  to  UUVs.  Doherty  [6] 
and  Prestero  [7]  both  employed  various  analytical  techniques  based  on  assumptions, 
geometry,  and  empirical  formulae  to  determine  the  hydrodynamic  coefficients  of  the 
Hydroid  REMUS100  UUV.  Prestero’s  work  dealt  with  a  typical  REMUS  while  Doherty 
derived  coefficients  for  the  “long-body”  REMUS  equipped  with  cross  body  thrusters. 
Both  authors  used  experimental  data  to  confirm  their  analytical  findings.  Chen  et  al.  [5] 
used  a  projective  mapping  method  to  estimate  the  dynamic  parameters  of  a  hovering  class 
UUV  similar  to  the  work  presented  here.  A  vision-based  system  is  used  to  capture 
several  images  of  a  UUV  during  specialized  maneuvers  to  track  the  vehicle  position.  Eng 
et  al.  [9]  used  a  free  decay  pendulum  technique  to  measure  the  hydrodynamic  coefficients 
directly. 

This  technique  differs  from  the  free-motion  techniques  mentioned  above  since  the 
vehicle  is  mounted  in  a  rig  that  constrains  the  vehicle  motion.  The  method  is  based  on  the 
classical  decay  test  with  damped  pendulum  motion.  A  UUV  is  mounted  to  the  pendulum 
rig  and  allowed  to  swing  freely.  It’s  planar  and  angular  motion  is  captured  on  video 
which  was  then  processed  to  determine  Cartesian  position  coordinates,  as  well  angular 
rate.  This  allowed  for  analytical  computation  of  hydrodynamic  parameters  using  the 
pendulum  dynamics,  and  not  the  free-motion  vehicle  dynamics,  limiting  applicability  to 
the  problem  at  hand.  Bahrke  [19]  used  analytical  as  well  as  on-line  statistical  estimation 
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to  estimate  the  speed,  steering,  and  diving  parameters  of  the  equations  of  motion  for  the 
NPS  AUV  II  vehicle.  A  Kalman  filter  based  parameter  estimator  is  used  to  estimate 
parameters. 

Parametric  models  are  infonned  by  the  physical  characteristics  of  the  system 
derived  from  first  principles  (e.g.,  the  generic  hydrodynamic  model  developed  by  Fossen 
[1]).  The  generic  hydrodynamic  model  informs  the  basis  functions  required  to  model  the 
system  response,  and  the  parameters  define  the  relative  weighting  of  these  basis 
functions.  Alternatively,  generic  sets  of  basis  functions  can  be  used  to  model  the  system 
dynamics.  One  such  example  is  neural  network  representations,  where,  for  example, 
radial  basis  functions  are  used.  System  identification  techniques  can  similarly  be  used  to 
learn  the  relative  weights  of  these  basis  functions.  These  techniques  are  referred  to  as 
non-parametric  approaches,  a  misnomer  since  the  model  is  in  fact  parameterized,  but 
these  parameters  do  not  have  physical  meaning.  Non-parametric  model  representations 
for  dynamical  systems,  such  as  neural  networks,  have  been  investigated  (e.g.,  Scarselli 
and  Tsoi  [10]).  According  to  Scarselli,  a  neural  network  with  sufficient  number  of  nodes 
is  capable  of  learning  and  representing  any  dynamic  system.  These  non-parametric 
models  are  less  desirable  since  they  lack  physical  meaning,  but  accurately  capture 
complex  input-output  relations  and  thus  find  application  is  simulator  development.  A 
neural  network  approach  has  been  applied  to  a  UUV  by  Juan  et  al.  [14].  The  motion 
model  of  a  Beaver  underwater  vehicle  was  developed  using  a  wavelet  neural  network 
while  the  thruster  model  was  developed  using  an  improved  radial  basis  function  neural 
network. 

C.  SCOPE  OF  THIS  WORK 

This  work  develops  a  non-linear  model,  based  on  the  rigid  body  dynamics  and 
hydrodynamic  force  and  moment  analysis  presented  by  Fossen,  for  a  hovering-class  UUV 
and  learns  the  model  parameters  using  various  system  identification  techniques.  This 
specialized  model  is  first  implemented  in  simulation  using  known  coefficients  (obtained 
from  Chen  et  al.  [5]),  which  is  used  to  verify  the  system  identification  and  online 
parameter  estimation  techniques  (e.g.,  error  convergence  and  parameter  estimate 
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convergence  to  the  true  values,  persistence  of  excitation,  etc.).  As  a  part  of  this  analysis, 
the  sequence  of  excitations  is  identified  to  learn  the  relevant  model  parameters.  Next, 
these  system  identification  techniques  are  applied  to  the  physical  system,  the  SeaBotix 
vLBV300,  to  obtain  the  parameters  for  the  simplified  hydrodynamic  model.  Position  and 
orientation  measurements,  obtained  from  an  instrumented  dive  tank  environment  and 
external  motion  capture  system,  are  used  to  learn  the  dynamic  model  of  the  UUV  in  real¬ 
time.  The  predicted  system  response  from  the  model  is  compared  to  the  true  response  of 
the  platform  to  validate  the  model.  Similarly,  a  probabilistic  estimation  technique  is 
evaluated  for  the  system  at  hand.  Finally,  a  non-parametric  model  is  developed  to 
evaluate  applicability  and  limitations  of  this  approach.  The  thesis  is  structured  as 
follows: 

Chapter  II  presents  a  brief  history  of  the  platform  of  interest,  the  SeaBotix 
vLBV300.  An  overview  of  its  current  and  potential  applications  is  presented  as  are 
appropriate  vehicle  specifications  provided  by  the  manufacturer. 

Chapter  III  presents  the  development  and  specialization  of  the  hydrodynamic 
(parametric)  model  for  a  hovering  class  UUV.  The  appropriate  assumptions  about  the 
UUV  and  its  assumed  operating  environment  are  presented  as  is  the  proposed  reference 
frames.  These  assumptions  allow  for  the  simplification  of  the  hydrodynamic  model. 
Next,  the  appropriate  and  complete  kinematic  and  dynamic  equations  of  motion  are 
developed  and  presented.  Finally,  a  propulsion  model  is  developed  for  the  individual 
thrusters  as  well  as  the  vectored  thruster  system  and  presented. 

Chapter  IV  discusses  the  various  system  identification  techniques  investigated. 
Recursive  methods  such  as  prediction  error  based  techniques  using  gradient  estimators  or 
least-squared  estimators  are  investigated  as  are  a  statistical  approach  based  on  Bayesian 
Inferencing,  and  a  non-parametric  approach  based  on  neural  networks. 

Chapter  V  presents  the  results  of  a  simulated  parameter  estimation  experiment 
based  on  known  parameters  from  previous  research  [5].  The  framework  for  the 
parametric  model  derived  in  Chapter  III  and  the  system  identification  techniques  derived 
in  Chapter  IV  are  validated  and  results  presented.  Also  presented  are  the  parameter 
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estimation  results  and  a  complete  hydrodynamic  model  for  the  vLBV300  vehicle 
developed  from  real  world  implementation  of  the  techniques  discussed  in  Chapters  III 
and  IV.  This  chapter  also  presents  a  verification  of  the  results  as  well  as  sequential 
procedure  for  detennining  the  dynamic  coefficients  of  a  UUV  using  on-line  parameter 
estimation  based  on  conclusions  obtained  from  the  experimental  data 

Chapter  VI  presents  conclusions  based  on  the  results  and  limitations  identified  by 
this  work.  It  also  presents  areas  for  further  refinement  or  research  and  conclusions  on 
this  work. 
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II. 


SEABOTIX  VLBV300  REMOTELY  OPERATED  VEHICLE 


A.  OVERVIEW 

The  vLBV300,  manufactured  by  SeaBotix,  Inc.  of  San  Diego,  California,  is  a 
MiniROV/Hovering  class  remotely  operated  vehicle.  It  is  operated  through  a  tether  and 
can  be  controlled  with  a  joystick.  Data  and  power  are  transmitted  through  the  tethered 
interface.  For  propulsion  and  maneuvering  the  vLBV300  uses  six  brushless  DC  thrusters, 
four  of  which  are  vectored  in  the  horizontal  plane  (surge,  sway,  yaw),  the  angles  of  which 
can  be  manually  adjusted,  while  the  remaining  two  are  fixed  in  the  vertical  plane  (heave, 
roll).  As  a  result,  the  vehicle  allows  control  in  5  degrees  of  freedom  (DOF)  when 
individual  thrusters  are  commanded,  or  4  DOF  when  controlled  via  the  joystick.  The 
vehicle  can  be  equipped  with  a  variety  of  sensors  including  a  controllable  FID  camera, 
rear  and  side  cameras,  sonars,  as  well  as  a  grabber  arm  for  in-water  intervention.  A 
typical  configuration  is  shown  in  Figure  1. 


Figure  1.  SeaBotix  vLBV300  miniROV  platform  (From  [15])  has  been  modified  to 

allow  tethered,  autonomous  operations. 

Based  on  its  specifications,  SeaBotix  recommends  the  vLBV300  for  a  variety  of 
applications  including  work  on  offshore  oil  or  gas  platforms,  coastal  and  inshore  surveys, 
maritime  security  of  ports,  harbors  and  vessels,  or  long  line  penetration  or  pipe 
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inspections.  The  vehicle  is  operated  in  tele-operated  mode  for  these  applications:  a 
human  pilot  controls  the  vehicle  via  the  joystick  in  the  surge,  sway,  heave,  and  yaw 
directions. 

In  the  current  research,  a  different  application  is  investigated:  autonomous  close- 
proximity  operations  in  the  presence  of  static  and  dynamic  obstacles  (i.e.,  human  divers). 
The  vLBV300  is  particularly  well  suited  to  proximal  operations  for  the  following 
reasons: 

•  The  vehicle  has  5  degrees  of  freedom,  including  surge,  sway,  heave,  yaw, 
and  roll.  Controllability  in  the  sway  direction  is  of  particular  importance 
for  proximal  operations  to  ensure  diver  safety. 

•  The  vehicle  is  light-weight  (i.e.,  easy  to  deploy  and  recover),  neutrally 
buoyant  (to  ensure  safe  operations  among  divers),  and  powerful  enough  to 
maneuver  among  divers  and  carry  considerable  payloads. 

•  The  vehicle  has  an  open-frame  architecture,  allowing  various  payloads  to 
be  integrated  with  relative  ease. 

•  A  computer  control  interface  has  been  developed  in  addition  to  the 
joystick  interface  that  allows  autonomous  operation  of  the  vehicle.  A 
high-level  control  interface  allows  computer  control  via  the  joystick 
commands,  while  a  low-level  control  interface  allows  individual  thrusters 
to  be  commanded. 

B.  VEHICLE  SPECIFICATIONS 

Specifications  for  the  vehicle  are  provided  by  the  manufacturer.  Table  1  lists 
general  specification  of  the  vLBV300  and  Table  2  lists  tether  specifications.  The  specific 


vehicle  setup  for  the  scope  of  this  work  is  discussed  in  the  next  section. 


Depth  Rating 

304.8  m 

Length 

0.625  m 

Width 

0.39  m 

Height 

0.39  m 

Diagonal 

0.55  m 

Weight  in  air 

18.09  kg 

Table  1.  General  characteristics  (From  [15]) 
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Diameter 

0.889  cm 

Length 

250  m  (nominal) 

Working  load 

978.6  N 

Breaking  strength 

6.86  kN 

Buoyancy 

neutral 

Table  2.  Tether  characteristics  (From  [15]) 


C.  CAVR  VEHICLE  CONFIGURATION 

The  vLBV300  is  a  miniROV  remotely  operated  vehicle  which  can  be  commanded 
via  a  high-  or  low-level  computer  interface,  resulting  in  a  tethered,  hovering-class 
autonomous  underwater  system  (THAUS).  As  configured  for  this  research,  the  vehicle 
weighs  20.9  kg  in  air;  is  0.625  m  long,  0.39  m  wide,  and  0.39  m  tall;  and  has  a  BlueView 
sonar  attached  in  addition  to  the  standard  payload  (tilt-controlled  forward-looking  camera 
and  LED  arrays).  This  payload  will  be  expanded  with  an  inertial  navigation  system 
(INS),  Doppler  velocity  log  (DVL)  and  GPS  unit  in  the  near  future. 

The  aft  pair  of  horizontal  thrusters  is  vectored  to  an  angle  of  45  degrees  from 
centerline  while  the  forward  pair  of  horizontal  thrusters  is  vectored  to  35  degrees  from 
centerline.  The  vertical  thrusters  are  angled  at  18  degrees  (see  Figure  2). 


A  right  handed,  body-fixed  frame  is  used  with  x  pointing  forward,  y  pointing 
right,  and  z  pointing  down  (see  Figure  (3)). 


Heave 

Figure  3.  Body  reference  frame  of  the  vLBV300 


D.  PROPULSION  MODEL 

As  discussed  in  Section  B,  the  SeaBotix  ROV  can  be  controlled  either  via  a  high- 
level,  joystick  interface  that  allows  for  four  DOF  control  or  via  low-level  pulse  width 
modulation  (PWM)  thruster  commands  directly  to  each  thruster  via  a  computer  program 
which  provides  for  five  DOF  control.  In  order  to  ensure  maximum  control  of  the  vehicle, 
for  this  research  only  the  low-level  commands  are  used. 

The  low-level  PWM  commands  range  from  -102  to  102  and  is  related  to  the 
thruster  RPM.  For  the  aft  thruster-pair,  the  positive  direction  is  defined  as  producing 
clockwise  rotation  of  the  propeller,  which  results  in  motion  in  the  positive  (ahead)  surge 
direction.  For  the  forward  thrusters,  a  positive  command  results  in  counter-clockwise 
rotation  and  therefore  also  thrust  in  the  positive  surge  direction.  A  positive  sway 
command  results  in  the  leading  edge  propellers  (right  hand  set  for  motion  to  the  right) 
spinning  in  the  positive  direction  while  the  trailing  set  spins  in  the  negative  direction.  A 
positive  command  to  the  vertical  thrusters  results  in  clockwise  rotation  and  motion  in  the 
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positive  heave  (downward)  direction.  Using  these  definitions  and  the  thruster  geometry, 
the  vehicle’s  propulsion  forces  are  defined  below.  For  ease  of  notation  the  thrust 
generated  by  the  horizontal,  forward,  starboard  propeller  is  referred  to  here  as  FS.  The 
horizontal,  forward,  port  propeller  thrust  is  referred  to  as  FP.  The  horizontal,  aft, 
starboard  propeller  thrust  is  referred  to  as  AS  and  the  horizontal,  aft,  port  propeller  thrust 
is  referred  to  as  AP.  The  thrust  generated  by  the  vertical  starboard  and  vertical  port 
propellers  are  referred  to  as  VS  and  VP,  respectively.  The  vectoring  angle  for  the  forward 
pair  of  thrusters  is  called  a  , ,  the  vectoring  angle  for  the  aft  pair  of  thrusters  is  called  ar , 

and  the  vectoring  angle  for  the  vertical  pair  of  thrusters  is  called  /3  (see  Figure  2).  Xprop 
is  the  total  force  generated  by  all  thrusters  in  the  surge  direction,  Yprop  is  the  total  force 
generated  by  all  thrusters  in  the  sway  direction,  and  Zprop  is  the  total  force  generated  by 
all  thrusters  in  the  heave  direction.  Kprop  is  the  total  moment  generated  around  the  body 
x-axis,  Mprop  is  the  total  moment  generated  around  the  body  y-axis,  and  Nprop  is  the  total 
moment  generated  around  the  body  z-axis. 


prop 


{^AP  *car  -  FP*c  oij  j  + 
AS  *  c  ar  -  FS  *  c  a  f  j 


(4) 


Yprop  ~  *  s  a  ^  +  AS  *  s  ar  j  - 


(fP  *  s  a  f  +  AP  *  s  ar  j 

(5) 

-(VL-VR)*sj3 

Zprop  =(VL  +  VR)*c/3 

(6) 

Kprop  =(VL  +  VR)*sP*LV 

(7) 

M  prop  ~  0 

(8) 
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N prop  ~  [^AP  *  c  ar  -  FP  *  c  a  ^  j  *  L//|  + 

( AS  *  c  ar  -  FS  *  c  a  , ■  J  *  LFl ,  + 

V  7  (9) 

*  s  ar  -  FP  *  s  ay  j  *  + 

|t.S  *  s  ar  -  FS  *  s  «y  j  *  LFl 2  + 


Note:  For  ease  of  notation  s*  =  sin(*),  c*=cos(*),  and  t*=tan(*).  ZF,  LHi,  and 
LH2  are  moment  arms.  LV  is  defined  as  the  distance  between  the  centroid  of  the  two 
vertical  propellers  while  LHi  is  the  distance  between  the  centroid  of  the  propellers  of  the 
aft  or  forward  pair,  and  LH2  is  defined  as  the  distance  between  the  centroid  of  the 
propellers  of  the  right  or  left  side  pair.  Mprop  is  zero  in  all  cases  since  pitch  cannot  be 
controlled  on  any  vehicle  considered  in  this  research. 

The  relationship  between  low-level  PWM  commands  and  a  measured  force  for 
individual  thrusters  is  determined  through  in -tank  testing  using  the  aft  thruster  pair  only 
and  a  Futek  USB210  strain  gage.  The  results  for  both  the  positive  and  negative  direction 
are  virtually  identical  for  commands  from  10  to  60,  as  would  be  expected  from  the 
performance  of  a  fixed  pitch  propeller.  Due  to  a  limitation  with  the  experimental  setup, 
surface  effects  affected  results  above  a  PWM  command  of  60  (i.e.,  the  60-102  range)  as 
can  be  seen  in  Figure  4.  The  vehicle  was  secured  to  an  anchor  point  outside  the  tank  via 
a  line  (with  the  strain  gage  integrated  into  this  line)  in  order  to  keep  the  strain  gage  dry. 
Because  the  line  rested  on  the  lip  of  the  tank  which  is  necessarily  several  inches  higher 
than  the  water  line,  the  vehicle  assumed  a  negative  pitch  (downward)  angle  during  every 
test.  This  resulted  in  the  aft  pair  of  thrusters  being  up  to  several  inches  closer  to  the 
surface  at  higher  PWM  commands  than  the  forward  thrusters.  This  difference  in  depth, 
and  therefore  water  pressure,  was  significant  enough  above  the  PWM  command  of  60  for 
the  aft  pair  of  thrusters  to  experience  significant  cavitation,  reducing  their  effectiveness. 
The  aft  pair  of  thrusters  was  used  for  the  positive  direction  testing  and  the  forward  pair  of 
thrusters  was  used  for  the  negative  direction  testing. 


14 


Force  Generated  by  a  Single  Thruster 


Figure  4.  Thruster  testing  showing  surface  effects 


Because  of  this,  a  regression  was  only  fit  to  data  for  PWM  commands  between  10  and  60. 
Since  the  planned  vehicle  operations  will  be  at  low  thrust  (and  speed),  this  thruster  model 
is  sufficient.  The  truncated,  experimentally  obtained  data  was  then  used  to  fit  a  second- 
order  polynomial  to  relate  PWM  values  to  generated  thrust. 

Thrust  =  0.006736PWM2  -  0.03366PWM  +  0.0684  (10) 

The  results  of  the  regression  are  compared  with  the  experimental  data  in  Figure  5. 
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Force  Generated  by  a  Single  Thruster 


Figure  5.  Bounded  propulsion  model  for  a  single  thruster  on  the  vLBV300  compared 

with  experimental  data 
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III.  GENERIC  MOTION  MODEL 


A.  FULL  EQUATIONS  OF  MOTION 

The  non-linear  dynamic  equations  of  motion  for  a  six  degree  of  freedom  (DOF) 
submersed  body  are  presented  by  Fossen  [1]  as 

Mu  +  C(v)v  +  D(v)u  +  g(rj)  =  r  /in 


The  vector  n  contains  the  body  velocities  and  angular  rates, 

u  =  [u,  v,  w,  p,  q, ) -f 


(12) 


where  u  is  velocity  in  the  surge  direction,  v  is  velocity  in  the  sway  direction,  w  is  velocity 
in  the  heave  direction,  p  is  the  angular  rate  about  the  x-axis  of  the  body  frame,  q  is  the 
angular  rate  about  the  v-axis  of  the  body  frame,  and  r  is  the  angular  rate  about  the  z-axis 
of  the  body.  All  velocities  in  this  work  are  presented  in  units  of  meters  per  second  and  all 
angular  rates  in  radians  per  second. 


The  vehicle’s  pose  (position  and  orientation)  is  described  by  the  vector  ij : 

n  =  [x,y,z,</>,0,yf  (13) 


where  x,  y,  and  z  are  the  vehicle’s  position  with  respect  to  an  inertial  reference  frame,  and 
<1 i>,  0,  y  are  Euler  angles  determined  using  the  zvx-convention.  All  positions  in  this  work 
are  presented  in  units  of  meters  and  all  Euler  angles  in  units  of  radians. 

Vehicle  pose  dynamics  are  related  to  body  velocities,  angular  rates,  and  the  pose 
itself,  by 

r/  =  J(q)u  (14) 


J  = 


J\ 


3x3 


(15) 


where 


A  = 


zy  cO 
sy  cO 
-s0 


-  s  y  c  <p  +  c  y  s  6  s  (j> 
c  y  c  (f>  +  s  <f>  s  9  s  y 
cd&(j) 


sy  s<t>  +  cy  c<j)s9 
-cys(j)  +  s9syc<j> 
c9c</> 


(16) 
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(17) 


J2 


1  sOtO  c<j)td 

0  c<j)  -sift 

0  scf)  led  of)  /  cB 


Matrices  M  and  C(v)  are  each  composed  of  two  parts:  a  contribution  from  the  rigid- 
body  dynamics  and  a  contribution  from  the  added,  or  apparent,  mass  and  inertia.  The 
rigid  body  dynamics  are  derived  by  Fossen  [1]  based  on  Newton’s  second  law  and 
Euler’s  first  and  second  axioms  to  that  law  while  the  added  mass  and  inertia  terms  are 
derived  from  an  energy  based  approach  resulting  from  Kirchoff  s  equations  as  discussed 

m[l]- 


The  inertia  matrix,  M,  is  the  sum  of  the  rigid  body  mass,  Mrb  ,  and  the  added  inertia 
matrix,  MA, 


M 

=  Mrb 

+  ma 
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(19) 


(20) 


where  m  is  the  vehicle  mass  in  air  and  xg,  y„,  and  zg  are  the  Cartesian  coordinates  in  the 
body  frame  of  the  vehicle’s  center  of  gravity.  Ix  is  the  vehicle’s  moment  of  inertia 
around  the  x-axis  of  the  body,  />  is  the  vehicle’s  moment  of  inertia  around  the  y-axis  of 
the  body,  and  Iz  is  the  vehicle’s  moment  of  inertia  around  the  z-axis  of  the  body.  These 
terms  represent  the  moment  of  inertia  due  to  un-coupled,  single  DOF  rotation  about  the 
axis  of  the  DOF  movement.  The  inertia  terms  with  more  than  one  letter  in  the  subscript 
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denote  the  vehicle’s  moment  of  inertia  due  to  coupled  motion.  For  example  Ixy  denotes 
the  vehicle’s  moment  of  inertia  around  the  vehicle’s  x-axis  due  to  rotation  around  the  y- 
axis. 


Fossen  defines  the  added  mass  coefficients  shown  in  Equation  (20)  as  describing 
the  pressure-induced  forces  and  moments  due  to  a  forced  harmonic  motion  of  the  body 
which  are  proportional  to  the  acceleration  of  the  body  [1].  For  example  Xv  is  the  force 
along  the  x-axis  due  to  an  acceleration,  v ,  in  the  y-direction.  Mathematically,  this  is 
expressed  as 


(21) 


C(v)  is  the  matrix  of  Coriolis  and  centripetal  terms  and  is  also  the  sum  of  a  rigid  body, 
Crb,  and  added  mass,  Ca,  components: 

C(v)±CRB{v)  +  CA(v)  (22) 


where 


Cm(p)  = 


3x3 

c. 


C 

c3 


(23) 


and 


C  = 


m(ygq  +  zGr ) 
~m(ygp  +  w) 
-m(zGp-v) 


-m(xgq  -  w) 
m(zGr  +  xgp ) 
-m(zGq  +u) 


-m{x  r  +  v) 
~ m(ygr-u ) 
m(xgp  +  ygq) 


(24) 


C2  = 


-m{ygq  +  zGr )  m(ygp  +  w)  m(zGp  -  v) 

m(x  q  —  w)  -m(zGr  +  xgp)  m(zGq  +  u) 


m{x  r  +  v) 


m(ygr-u ) 


-m(xgp  +  ygq) 


C3  = 


0  ~^rz<l~  1 xzP  ^zr  Irzr  + 1 xyP  ~ lyQ 

lyztf  ^xzP  ~  ^zr  ®  xzr  ~  I xiP  ^  I xP 

—Iyzr  ~  IxrP  +  I xzr  +  ^xrQ  ~  ^xP  0 


(25) 


(26) 
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0  0  0  0  -a3  a2 

0  0  0  a3  0  —ax 

0  0  0  -a2  ax  0 

0  -a3  a2  0  — Z>3  b2 

a3  0  -a,  b3  0  — 

-a2  Oj  0  -b2  bx  0 

where 

«1  =  XuU  +  XvV  +  xww  +  XPP  +  xfl  +  'V 
a2  =  X,u  +  7.v  +  Kw  +  Ypp  +  Yq  +  Lr 

«3  =  +  V  +  Z«,W  +  +  Z,//  +  Z,'' 

hi  =  x.  u  +  Y.V  +  Zpw  +  KpP  +  K.q  +  K,r 
b2  =  X.u  +Y.v  +  Z.w  +  K.p  +  M.q  +  Mr 
b3  =  X.u  +  Y.v  +  Z.w+  K.p  +  M.q  +  N.r 

The  damping  matrix,  D(v),  is  the  sum  of  radiation-induced  potential  damping  due 
to  forced  body  oscillation.  Dp,  linear  skin  friction  due  to  laminar  boundary  flow  and 
quadratic  skin  friction  due  to  turbulent  flow,  Ds,  wave  drift  damping,  Dw,  and  damping 
due  to  vortex  shedding,  DM, 

D(v)±Dp(v)  +  Ds(v)  +  Dw  (u)  +  Dm  (u)  (29) 

Radiation  induced  damping,  DP  is  often  referred  to  as  potential  damping  and  is  a 
result  of  a  body  being  forced  to  oscillate  with  the  excitation  frequency  of  waves 
encountered  which  results  in  added  mass,  damping,  and  restoring  forces.  While  generally 
negligible  compared  to  other  forces  at  great  depth  for  underwater  vehicles,  it  is  of  more 
concern  for  surface  vehicles. 

Damping  due  to  skin  friction  is  a  function  of  the  vehicle’s  exterior  make  up  as  well  as 
speed.  Low-speed,  laminar  flow  results  in  a  low  frequency  contribution  while  turbulent 
flow  results  in  a  high  frequency  contribution. 

Wave  drift  damping  is  only  of  significance  for  surface  vessels  advancing  into  waves. 
Damping  due  to  vortex  shedding  occurs  in  a  non-viscous  fluid  and  is  a  function  of  the 
speed  at  which  the  vehicle  moves,  the  density  of  the  water  in  which  it  is  operating,  the 
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projected  cross  sectional  area  and  the  vehicle’s  Reynolds  number.  This  is  the  most 
significant  of  the  damping  contributors  for  an  underwater  vehicle.  The  general  equation 
for  this  damping  is  further  explored  by  Equation  (45).  g(q)  is  the  vector  of  restorative 
forces  and  moments, 


g(> 7)  = 


(W-B)s0 

-(W-B)c0s<l> 

-{W-B)c0c</> 

~(yGW - yBB) c0c</>  +  (zGW  -zbB)c0s</> 
(zGW  -  zbB) s0  +  (xgW  - xbB ) c0c<!> 
-(xGW-xBB)c0s0-(yGW-yBB)s0 


(30) 


where  W  is  the  weight  (in  air)  of  the  vehicle,  W=  mg,  B  is  buoyant  force  produced  by  the 
fully  submerged  vehicle,  and  xg,  yg,  and  zg  are  the  Cartesian  coordinates  in  the  body 
frame  of  the  vehicle’s  center  of  gravity  while  Xh,  y/>,  and  z/,  are  the  Cartesian  coordinates 
in  the  body  frame  of  the  vehicle’s  center  of  buoyancy.  All  gravitational  forces  and 
moments  are  assumed  to  act  through  and  around  the  vehicle’s  center  of  gravity  while  all 
buoyancy  forces  and  moments  are  assumed  to  act  through  and  around  the  vehicle’s  center 
of buoyancy. 


The  vector  of  the  forces  and  moments  exerted  on  the  vehicle  by  the  thrusters  is  r. 
Recall  from  Equations  (4)  -  (9)  that  Xprop,  Yprop,  and  Zprop  are  forces  in  the  translational 
DOF  while  Kprop,  Mprop,  and  Nprop  are  moments  in  the  rotational  DOF. 

r  =  [Xprop  Yprop  ^prop  K prop  ^  prop  prop  ]  (31) 


The  sign  convention  for  the  forces  that  make  up  r  follow  the  right  handed  reference 
frame  sign  convention  established  for  the  vehicle  in  Chapter  II.  That  is,  Xprop  is  positive 
when  it  generates  a  force  in  the  positive  surge  direction  (forward).  The  force  Yprop  is 
positive  when  it  generates  a  force  in  the  positive  sway  direction  (right).  The  force  Zprop  is 
positive  when  it  generates  a  force  in  the  positive  heave  direction  (downward).  The 
moment  Kprop  is  positive  when  it  generates  a  right-handed  moment  around  the  x-axis. 
The  moment  Mprop  is  positive  when  it  generates  a  right-handed  moment  around  the  y- 
axis.  The  moment  Nprop  is  positive  when  it  generates  a  right  handed  moment  around  the 
z-axis. 
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B.  ASSUMPTIONS  AND  SIMPLIFIED  EQUATIONS 


In  order  to  model  the  non-linear  dynamics  of  the  UUV,  Equations  (11)  and  (14) 
must  be  solved  for  v  and  ;/.  Several  assumptions  are  made  about  the  characteristics  of  the 
UUV  to  simplify  the  components  of  Equation(l  1).  The  vehicle  is  assumed  to  be 
neutrally  buoyant  (B=W)  and  the  tether  is  assumed  to  have  no  effect  on  vehicle  motion  or 
trim.  The  x  and  y  coordinates  of  the  center  of  gravity  and  center  of  buoyancy  are 
assumed  to  coincide.  Also,  since  the  specific  mass  distribution  of  the  vehicle  is  not 
known,  a  uniform  distribution  of  mass  is  assumed  [5,  6]  and  the  components  of  the  inertia 
tensor  that  are  off  the  main  diagonal  are  neglected.  The  moments  of  inertia  are  estimated 
using  the  equation  for  moment  of  inertia  of  a  solid,  rectangular  cuboid 


Ix  =  —m(height2  +  width1) 


IY  = — m(  height1  +  length2) 


Iz  =  —m(  length^  +  width2) 


(32) 

(33) 

(34) 


and  result  in  Ix=  0.559  kgm2,  IY  =  0.953  kgm2,  and  Iz  =  0.90  kgm2  for  the  vLBV300. 
These  assumptions  further  simplify  the  equations  of  motion  to  yield: 


gin)  = 


o 

o 

o 

~B(zg  -zB)cOs(j) 
B(zg-zb)s0 
0 


(35) 


Fbk  =  B(zg  -ZB)c0St/> 

(36) 

Fbm  =  B(ZG  ~zb) 

(37) 
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m 


where 


000  mzG  0 
0  m  0  —  mzG  0  0 

0  0  m  0  0  0 

0  ~mzG  0  I x  0  0 

mza  000  I Y  0 

o  ooo  o  /z 


Cm(p) 


03X3  Cl 

c2  c3 


C2  = 


mzGr 
—mw 

m(zGp-v) 

-mzGr  mw  m(zGp  —  v) 
-mw  —mzGr  m(zGq  +  u ) 
mv  -mu  0 


mw  —mv 

mzGr  mu 

-m(zGq  +  u)  0 


Q  = 


0 

-v 


V 

o 


-IYq 

Ixp 


IYq  -Ixp  0 


(38) 


(39) 


(40) 


(41) 


(42) 


Furthermore,  the  vehicle  is  assumed  to  have  three  planes  of  symmetry,  will  be 
fully  submerged  for  all  motion,  and  will  only  move  at  low  speeds.  Therefore,  (see  [1] 
Equations  (2.17)  and  (2.23))  the  following  simplified  terms  are  obtained: 


M A  =  -diag 


X.,Y.,Z.,K.M.,N. 
u ’  v’  w’  p  q  r 


(43) 
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-Mqq 
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(44) 


Assuming  in  addition  that  the  vehicle  is  performing  only  uncoupled  motion  with 
complete  turbulent  flow,  the  damping  term  is  in  the  form: 
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(45) 


D(u)  diagj^  A;/|U|  |(/|,E(;|V|  |v|,ZW|W|  |w|,AT^|^|  |^r|, A^r|r|  |r|  j 

Combining  all  assumptions  yields  the  following  simplified  dynamic  equations  for 
non-coupled  motion  for  the  THAUS  UUV: 

Xprop  =  (/n  ~  Xu  )«  -  (m  -  YM  +  (m  ~  Z„)HV  ~  Xu\u\U  H 

Yprop  =  ( m  -  YJv  -  (m  -  Zw)wp  +  (m  -  X .  )ur  -  F|v|v  |v| 

Zprop  =  ( m  - ZJw - (m  - X. )uq  +  (m  -  Y.)vp  - Zwjwjw | w| 

Kprop  =  (4  ~ Kp)p  +  (Yi- ZJwv  +  (7Z -IY+N-  +M. )rq 

-KMp\p\-FBK  (46) 

M ProP  =  {I y-M.)q  +  (Z,,  -  X,  )uw  +  (/Y  -  4  +  N,  -  K.  )rp 
-MMq\q\  +  FBM 

Xprop=Vz-Ni)r  +  {Xi,  —  YY)uv  +  (4  ~IX  +  Kp -M.)rp 
- N  i  ,r\r\ 

The  validity  of  these  assumptions,  based  on  the  accuracy  of  the  learned  model  compared 
to  the  measured  response,  will  be  discussed  in  Chapter  VI. 
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IV.  SYSTEM  IDENTIFICATION 


A.  PARAMETER  ESTIMATION  FOR  A  STATIC  SYSTEM 

Consider  a  regression  model  of  the  form  [8] 

yi=ri6l+tf01+...+(p”G"'=tf0 


(47) 


where  yit  i=l,2,...,n,  is  a  series  of  observations,  (pt  =  \jp]  <pf  ...  $"]  are  the  known 
regressor  functions  (which  are  functions  of  the  states  and  controls),  and 
©  =  \d  8~  ...  8"’  are  the  unknown  parameters  (to  be  estimated).  The  goal  of  online 

system  identification  is  to  recursively  estimate  the  unknown  parameters  0  according  to 
some  adaptation  law  as  new  information  becomes  available. 

Several  methods  of  parameter  estimation  are  available  to  perfonn  the  necessary 
recursive  estimation,  including  the  gradient  estimator,  the  recursive  linear  least  squares 
(RLLS)  estimator,  a  statistical,  Bayesian  approach,  and  a  neural  network  based  approach. 


1.  Recursive  Linear  Least  Squares 

The  goal  of  the  linear  least  squares  estimator  is  to  minimize  the  square  of  the  error 
between  the  complete  history  of  the  measured  and  modeled  response  of  the  system  with 
respect  to  0,  as  in  [2]  and  [8]. 

1 

min  V (©,  n)  =  min  —  £  (  y.  -  <pj®)2  (48) 

0  0  2  i=l 


T 

To  accomplish  this,  let  Y  =  [y]  y2  ...  yn  ]  be  the  set  of  measurements  and  define  the 
error  terms  sj  =  yj  -  <p'  ® ,  where  ©  is  an  estimate  of  the  unknown  parameter  vector.  This 
can  be  written  in  vector  notation:  E  =  [f,  s2  ...  sn^  .  Finally,  let 


T 

<Pl 


(t  = 


WA 


(49) 


Then  the  loss  function  in  Equation  (48)  can  be  written  as: 
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(50) 


1  11  O  1  T  1  II  ||2 

V(d,  n)  =  -  s  e-  =  - . E tE  =  -  £ 

2  /=i  2  2"  11 


where  E  =  Y  -  O0  .  This  loss  function  is  minimized  by: 

0*  =  (ct7®)1  q/y  (51) 

which  requires  that  ®r®  is  non-singular  to  yield  a  unique  solution.  In  a  similar  manner, 
the  recursive  linear  least  square  (RLLS)  estimator  is  defined  as 

0«=0h  +  ^,-^h)  (52) 

where  the  filter  gain  is  given  by: 

Ki=Pi_l<pi(l  +  ^Pi_l<pi)~1  (53) 

k 

and  P  =(I -Kxp] )P  , .  To  ensure  a  unique  solution  to  RLLS  estimation  ®r®  =  £  w <p[ 

i= 1 

must  be  of  full  rank,  or  non-singular. 


2.  Neural  Network 

Computational  neural  networks  (NN)  are  modeled  after  the  understanding  of 
biological  nervous  systems  and  modeled  as  a  dense  web  of  elements  or  nodes  with 
weighted  connections  that  may  be  adapted,  or  trained,  during  use  of  the  network  to 
improve  performance  [16].  The  strength  of  the  NN  lies  in  its  ability  to  process  many 
computational  tasks  in  parallel  instead  of  in  a  serial  manner  as  done  by  most  computers, 
in  its  inherent  robustness,  and  in  its  ability  for  adaptation  or  learning.  The  large  number 
of  interconnected  nodes  allows  for  several  competing  hypotheses  to  be  investigated  and 
pursued  simultaneously  by  the  network.  Because  a  NN  is  composed  of  a  large  number  of 
interconnected  nodes,  damage  to  one,  or  even  several,  may  not  completely  destroy  the 
system  as  would  happen  in  a  serial  computer.  Neural  networks  using  variable  weights 
between  nodes  and  a  learning  strategy  also  have  the  ability  to  learn  and  improve 
performance  during  operation  or  after  damage.  This  ability  also  provides  for  a  degree  of 
robustness  in  the  presence  of  variability  of  the  processing  element. 
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Neural  networks  are  of  interest  to  this  work  because  of  this  ability  to  learn  and 
adapt  online  as  well  as  their  ability  to  accurately  model  complicated  non-linear  systems. 
While  the  end  results  is  a  non-parametric  “black  box”  input-output  model  that  does  not 
afford  the  same  physical  understanding  of  a  system  as  the  hydrodynamic  model  by 
Fossen  does,  the  NN  model  can  provide  a  highly  accurate  simulation  that  is  able  to 
account  for  the  nonlinearities  inherent  to  modeling  an  underwater  vehicle. 

The  simplest  type  of  neural  network  sums  N  inputs  of  x,  which  commonly 
represents  a  vehicle’s  state  or  control  input  such  as  thruster  commands,  x,  is  multiplied 
by  its  respective  weights,  w-b  and  the  network  passes  the  results  through  a  non-linear 
function,  /  with  a  delay  described  by  0  as  presented  by  Lippmann  [16].  The  non-linear 
function  is  called  the  “basis”  function  and  helps  define  the  architecture  of  the  model  in 
the  same  manner  as  the  regressor  does  in  a  parametric  representation.  A  visual 
representation  is  given  in  Figure  6. 

y  =  f(JjWix-e)  (54) 

i=0 


Nodes 


Figure  6.  A  basic  neural  network  showing  the  network  output  expressed  as  a  sum  of 

weighted  inputs. 

The  goal  of  a  large  group  of  nodes  becomes  to  determine  which  of  a  large  set  of 


functions  is  best  representative  of  an  unknown  input  function  or  series  of  input  functions, 
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as  well  as  non-parametric  uncertainties.  Non-parametric  uncertainties  typically 
encompass  measurement  noise,  equipment  failure,  enviromnental  disturbances  or  other 
uncertainties  [21]. 

In  system  identification  terms,  as  presented  in  Section  VI.A,  this  is  expressed  as: 

y.=0T(pi+si  (55) 

where  s.  represents  the  approximation  error  of  the  non-parametric  uncertainties.  In  order 

to  identify  the  parameters  contained  in  61  in  the  face  of  these  uncertainties,  a  set  of  basis 
functions,  (pt ,  which  in  the  RLLS  context  were  informed  by  the  parametric  model  itself, 

must  be  chosen.  These  basis  functions  are  chosen  such  that  the  error  becomes  small 
over  the  expected  area  of  operation  of  the  system  being  identified. 

Physically,  this  means  that  a  series  of  basis  functions  are  tested,  along  with 
weights  for  each  portion  of  the  basis  function  until  the  basis  function  set  and  associated 
weights  are  found  that  most  closely  represents  the  uncertain  input-output  relation. 
Lippmann  presents  an  adaptive  classifier  along  these  lines  (see  Figure  7)  [16]  that  takes  a 
series  of  inputs  and  computes  a  series  of  performance  related  scores  that  reflect  the 
accuracy  of  each  basis  function  examined.  These  scores  are  passed  on  to  a  second  stage 
where  only  the  best  score  is  passed.  This  best  score  is  then  used  to  adapt  the  weights  of 
the  network. 
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Figure  7.  Diagram  of  an  adaptive  neural  network  (After  [16]) 


“Training,”  or  adaptation  of  the  network,  can  take  place  online  as  in  Figure  7,  or 
prior  to  using  the  network.  Input  data  is  fed  to  the  network  in  concert  with  known  output 
data.  The  network  then  trains  itself  according  to  various  methods  and  a  set  of  coefficients 
are  produced.  Two  of  the  most  well-known  adaptation  schemes  are  the  steepest  decent 
algorithm,  or  error  back  propagation  (EBP),  and  the  Gauss-Newton  algorithm.  The  EBP 
method  is  presented  in  [17]  as 

wM=Wi+agi  (56) 

where  wt  is  the  weight  assigned  to  a  particular  node  at  step  i,  a  is  the  step  size,  gi  is  the 
first-order  derivative  of  the  total  error  function  E(x,w),  and  wi+[  is  the  updated  node 

weighting.  Further  notation  in  [17]  includes:  p  is  the  index  of  training  patterns  from  1  to 
P,  m  is  the  index  of  outputs  from  1  to  M,  and  N  is  the  total  number  of  weights. 

The  Gauss-Newton  method  uses  the  Jacobian  matrix,  J,  to  relate  the  gradient  g  to 
a  vector  of  errors  e. 
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(58) 

(59) 


This  results  in  the  Gauss-Newton  learning  algorithm  to  recursively  update  the  weights  of 
the  network. 

wi+i  =  wi-{JUi)~ljiei  (6°) 


Where  the  EBP  is  stable  it  converges  slowly  and  while  the  Gauss-Newton  method 
converges  rapidly,  it  is  unstable.  The  training  method  used  in  this  work  is  the  Levenberg- 
Marquardt  method  which  is  a  combination  of  the  steepest  descent  algorithm  and  the 
Gauss-Newton  algorithm.  The  Levenberg-Marquardt  method  combines  the  strengths  and 
omits  the  weaknesses  of  the  EBP  and  Gauss-Newton  methodology.  It  converges  rapidly 
and  is  stable  [17].  A  combination  coefficient,  // ,  is  introduced  and  the  learning  method 
becomes 


w,+l  =  w,  -  (J'i  J i + a/)1  J,ei 


(61) 
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where  /  is  the  identity  matrix.  The  combination  coefficient  is  used  to  switch  between  the 
EBP  and  Gauss-Newton  methods  during  training.  As  ju  approaches  zero  and  the  Gauss- 
Newton  method  is  used.  When  it  is  large,  the  EBP  method  is  used. 

The  computer  program  MATLAB  provides  a  NN  toolbox  for  the  creation  and 
training  of  various  types  of  NN.  For  this  work,  a  nonlinear  input-output  autoregressive 
neural  network  (NARX)  is  created  and  trained  with  the  Levenberg-Marquardt  method. 
Ten  nodes  were  used  along  with  four  delays  in  order  to  maximize  network  performance 
while  balancing  the  required  computational  infrastructure  based  on  the  principle  that  the 
number  of  nodes  and  delays  determines  how  well  the  neural  network  can  approximate  the 
system  being  modeled  [10].  The  network  architecture  is  shown  in  Figure  8. 


Figure  8.  Network  architecture  of  a  NARX  created  in  Matlab 

Using  this  toolbox  and  the  NARX  network  architecture  ,  these  methods  can  be 
applied  to  system  identification  (with  the  intent  to  learn  complex  input-output 
relationships).  In  terms  of  the  regression  model  of  Section  IV. A  known  data,  (e.g,  the 
regressor  data  contained  in  <p ),  can  be  fed  in  along  with  known  output  data  y.  The 
network  will  then  find  an  appropriate  basis  function  and  train  the  coefficients  to  closely 
match  input  to  output.  Another  benefit  to  the  NARX  architecture  over  other  architectures 
(including  the  simple  non-linear  input  output  model,  or  NIO)  is  that  past  values  of  y  are 
fed  back  to  the  network  if  they  are  available  to  provide  more  accurate  error  estimation, 
online  training,  and  an  overall  better  estimate. 

The  major  disadvantage  of  a  NN  model  of  a  physical  system  is  that  the  resulting 

model  does  not  have  any  physical  meaning.  As  discussed  before,  the  end  result  of  NN 
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estimation  is  a  non-parametric  (aka  “black  box”)  representation.  The  basis  functions  and 
weights  assigned  to  them  do  not  provide  any  physical  insight  into  the  system  as  Fossen’s 
hydrodynamic  coefficients  do.  So,  while  the  resulting  model  may  describe  vehicle 
motion  extremely  well  it  cannot  be  used  to  infer  other  physical  characteristics  of  the 
vehicle  and  must  be  validated  by  direct  comparison  of  simulation  to  actual  results.  The 
value  of  the  NN  lies  in  its  ability  to  leam  complex,  non-linear  systems  and  account  for 
these  non-linearities  through  the  training  process  whereas  the  parametric  representation 
relies  on  assumptions  and  simplifications  that  do  not  always  hold  true.  In  the  end,  a  very 
accurate  model  of  a  complex  system  can  be  produced  with  disproportionately  small  time 
and  effort. 


3.  Gradient  Estimator 

The  gradient  estimator  (GE)  is  the  simplest  of  the  various  on-line  estimation  tools 
[2]  and  is  a  prediction-error-based  method.  The  prediction  error  is  defined  as 


=  y,  -  y, 


(62) 


©  =  0-0 


(63) 


where  y  is  the  predicted  output  of  Equation  (47)  at  time  step  i,  ©  is  a  vector  of  the 
estimated  parameters,  and  ©  is  the  parameter  estimation  error  .  The  gradient  estimator 
works  by  updating  the  parameters  contained  in  ©  in  the  opposite  direction  of  the  gradient 
of  Equation  (62),  the  instantaneous  prediction  error,  so  that  e,  is  reduced. 

d[eje,] 


©  =  “A) 
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(64) 


where  po  is  the  estimator  gain.  Equation  (64)  can  be  re-written  in  terms  of  ©  and  cp  as 

0  -  ~p{)(p]ei  (65) 


and  in  terms  of  the  parameter  estimation  error  0  as 


0  =  -p0(pT(pS 


(66) 


Slotine  [2]  proves  the  stability  of  Equation  (65)  using  the  Lyapunov  candidate  function 
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V  =  ©r0 


(67) 


and  gives  its  derivative  as 

V  =  — 2 pJdT (pT (f®  <  0  (68) 

which  implies  that  the  GE  is  stable  and  the  magnitude  of  the  squared  parameter  error  is 
always  decreasing.  To  guarantee  that  the  GE  will  converge  to  the  true  value,  however, 
the  signal  used  to  excite  it  must  be  sufficiently  exciting.  A  signal  is  said  to  be  sufficiently 
exciting  if  it  possesses  the  appropriate  level  of  persistence  of  excitation  (PE).  That  is,  the 
input  signal  must  be  sufficiently  rich  in  frequency  content  (normally  accomplished  by 
using  multiple  sinusoids  as  inputs)  to  excite  the  relevant  modes.  Simply  put,  a  signal 
composed  of  several  sinusoids  is  capable  of  estimating  more  parameters  than  a  signal  of 
one  sinusoid  (this  concept  is  explored  in  depth  in  Section  IV. C).  For  this  work,  PE  is  not 
a  limiting  condition  since  multiple  sinusoidal  input  signals  may  be  applied  to  the  vehicle 
to  excite  it. 

For  the  gradient  estimator  PE  is  also  essential  to  ensure  the  robustness  of  the 
estimator.  Robustness  is  the  ability  of  the  estimator  to  maintain  reasonably  good 
parameter  estimation  in  the  presence  of  parametric  uncertainties  such  as  parameter  time- 
variation,  measurement  noise,  and  disturbances.  If  the  input  signal  is  not  sufficiently  PE 
then  the  estimator  values  may  diverge  even  without  the  presence  of  noise  or  uncertainty 
[2]. 

The  quality  of  the  estimates  produced  by  the  gradient  estimator  also  depends  on 
the  rate  of  parameter  variation,  the  level  of  non-parametric  uncertainties,  and  the 
magnitude  of  the  estimator  gain  [2], 

The  rate  of  variation  of  the  parameters  to  be  estimated  affects  the  robustness  of 
the  estimator  as  well.  If  the  parameters  vary  rapidly  the  estimator  will  have  a  much 
harder  time  accurately  estimating  the  parameters  and  converging  to  the  true  value  [2]. 
The  hydrodynamic  parameters  to  be  estimated  by  this  work  are  not  expected  to  be  time- 
varying  so  this  limitation  of  the  gradient  estimator  does  not  come  into  play  here.1 

1  The  parameters  have  been  experimentally  observed  to  vary  with  velocity.  In  order  to  compensate  for 
this  some  time -variance  may  be  assumed  later.  This  assumption  will  be  revisited  in  Chapter  VI. 
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The  presence  of  noise  or  un-modeled  disturbances  in  the  vehicle  dynamics  and 
measurements  play  a  large  part  in  the  accuracy  of  the  results  produced  by  the  GE.  This  is 
illustrated  by  adding  a  disturbance  term,  <4  to  y,  in  Equation  (62)  resulting  in 

ei=yi-(yi  +  di)  (69) 

Substituting  Equation  (69)  into  Equation  (65)  shows  the  effect  of  the  disturbance 

®  =  -p0<pT[yl-(yl+di)]  (VO) 

which  in  turn  affects  parameter  error  through  Equation  (63).  That  is,  for  a  large 
disturbance,  the  parameter  estimation  error  becomes  large.  For  this  work,  which  involves 
noise  in  both  the  measurement  noise  from  the  motion  capture  system  and  process  noise 
from  the  highly  non-linear  nature  of  a  body  moving  in  a  non-viscous  fluid,  disturbances 
will  greatly  degrade  the  accuracy  of  a  gradient  estimator. 

The  effect  of  the  magnitude  of  the  estimator  gain  is  easily  seen  in  Equation  (64). 
A  larger  gain  results  in  a  faster  rate  of  convergence,  but  just  as  in  optimization,  increasing 
the  gain  or  step  size  only  improves  performance  up  to  a  certain  point.  Increasing  the  gain 
to  too  high  a  level  can  result  in  oscillation  and  very  slow  convergence  to  the  true 
parameters.  Since  the  estimator  gain  is  defined  by  the  user  this  factor  does  not  limit  the 
use  of  a  GE  for  this  research. 


4.  Bayesian  Filtering 


Bayes’  rule  provides  a  method  to  compute  a  posterior  probability  from  a  set  of 
given  prior  probabilities.  It  states  that  given  a  collection  of  k  mutually  exclusive  and 
exhaustive  events  A  i,  A?,  ...,Ak  that  have  their  own  prior  probability  P(Ai)(i=l,...,k),  and 
any  other  event  B ,  with  P(B)  >  0,  the  posterior  probability  of  Aj  given  B  has  occurred, 
P(Aj  |  B) ,  can  be  expressed  as  [15] 


P(AJ  |  B )  = 


P{B\Aj)P{Aj) 

fjP{B\Aj)-P{Aj) 

i=l 


j  =l,.-,k 


(71) 
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This  concept  can  easily  be  applied  to  parameter  estimation.  Thrun  et  al.  [3] 
present  an  adaptation  that  uses  measurement  and  control  input  data  to  represent  and 
evolve  the  belief  state:  the  probability  of  state  x;  given  all  the  available  information.  The 
belief  about  a  state,  such  as  hydrodynamic  parameters  being  estimated,  at  a  time  i  is 
calculated  from  the  belief  at  time  t- 1 ,  the  control  input  ut  at  time  t,  and  the  measurement, 
yu  For  example,  assuming  that  xt  has  a  Gaussian  distribution,  the  belief  state,  b,  can  be 
described  by  the  first  and  second  moments  of  the  distribution,  mean,  /j ,  and  variance, 

a2 . 

b(xi)  =  —f=e  (72) 

<jxJ2n 

Bayes  rule  provides  an  algorithm  to  recursively  calculate  the  new  belief  state 
given  the  dynamic  and  measurement  models,  as  well  as  the  previous  belief  state  and  the 
measurement  (e.g.,  [3]).  It  consists  of  two  steps.  The  first  step  is  called  the  prediction 
and  it  calculates  the  belief  for  state  xh  called  b(xt)  ,  based  on  the  belief  of  the  previous 
belief-state,  xt-i,  and  the  control  input  a,. 

i 

b(xt )  =  J p(xt  |  xi_l,ui)b(xi_l)dxi_l  (73) 

o 

The  second  step  is  referred  to  as  the  measurement  update.  Here  the  Bayes  filter 
uses  the  belief  calculated  by  the  prediction  step  in  Equation  (58)  and  multiplies  it  by  the 
probability  that  measurement  z,  was  observed.  Each  hypothetical  posterior  state  is  treated 
this  way  and  because  the  resulting  product  may  not  integrate  to  one  the  nonnalization 
factor  77  is  used  to  ensure  a  proper  probability  is  returned. 

b(x,)  =  ?ip(zl\xl)b(xl)  (74) 

In  order  to  use  Equations  (73)-(74)  recursively  b(xo)  must  be  initialized  at  i  =  0.  Thrun  et 
al.  suggest  that  if  xo  is  precisely  known  then  b(xo)  should  be  initialized  as  a  point  mass  on 
the  correct  value  of  xo  with  zero  probability  elsewhere.  If  nothing  is  known  about  xo, 
however,  a  uniform  distribution  should  be  used  for  b(xo)  [3].  Thrun  et  al.  then  gives  the 
recursive  form  of  Equations  (73)-(74)  as 
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b  (x, )  =  J  p(xi  |  xi_l,ut)b(xi_l)dxi_l 


(75) 


b(xi)  =  T]p(zi\xi)b(xi) 


(76) 


Practical  implementation  of  Equations  (75)-(76)  requires  the  initial  belief  b(x0),  the 
measurement  probability  p(z[  |  x,) ,  and  the  state  transition  probability  p(xt  \  m.,xm )  .  The 

choice  of  the  probability  distributions  for  these  required  probabilities,  as  well  as  the 
system  being  modeled,  shapes  the  type  of  Bayesian  filter.  For  example,  assuming  a 
linear  system,  an  initial  Gaussian  distribution,  that  all  three  probability  distribution 
functions  are  zero  mean,  Gaussian  white  noise  turns,  and  that  the  dynamic  and 
measurement  equations  are  linear  the  Bayesian  filter  into  a  Kalman  filter 

For  the  problem  at  hand,  the  system  has  non-linear  dynamics  described  by: 

*/+l  =f(xi,U„G>i) 

co,~N(0,Q )  1  } 


where  x/+1is  a  vector  of  states  at  time  i+1,  x,.  is  a  vector  containing  the  states  to  be 


estimated  at  time  step  i,  ui  is  the  control  input  at  time  step  i,  and  (ot  is  the  zero  mean, 


white  noise  with  covariance  Q  associated  with  the  process  at  time  step  i.  The  system  has 
measurements  zp 


z,  =  h(xt,vt) 
K~7V(0,7?) 


(78) 


where  vk  is  the  zero  mean,  white  noise  measurement  noise  with  covariance  R  at  time  step 

i.  One  implementation  of  Bayes  rule  for  non-linear  systems  is  to  linearize  the  system 
around  trim  points  and  to  apply  the  Kalman  filter  machinery,  resulting  in  different  types 
of  filters  examples  of  which  are  the  extended  Kalman  filter,  unscented  Kalman  filter,  and 
the  particle  filter  [3].  To  do  this  the  concept  of  the  Jacobian  is  used.  The  Jacobian  is 
defined  as  the  partial  derivative  of  a  process  with  respect  to  each  of  the  variables  that 
inform  that  process.  Now,  define  the  following  i  x  j  Jacobian  matrices 
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A  a! 

J  dx: 


dm, 


(4,«t,o) 


a  dh, 

H°-aP, 


v  s^i 

11  dv, 


(xt,Uk,0) 


(4.0) 


(4.0) 


(79) 


Remembering  that  all  noise  distributions  are  assumed  to  be  zero  mean,  Gaussian 
white  noise  and  using  an  initial  state,  xq,  and  an  initial  covariance  matrix,  Pq,  the 
prediction  step  is: 


(80) 


where  jc!lM  is  the  predicted  state  at  time  i  given  the  state  information  from  the  previous 
step  and  /)|M  is  the  predicted  covariance  at  time  i  based  on  the  covariance  information 
from  the  previous  step.  Next  is  the  measurement  step 

yt=zi-h(xpf_  i)  (81) 


which  consists  of  the  measurement  z;  and  gives  as  its  output  the  updated  measurement 
error  y. .  The  Kalman  gain,  A),  is  then  calculated  as 


+  VtRV?Tl 

and  the  estimate  xt  is  updated  by 

%=xAi-x+Ki(zi-KxiV_l,Q)) 

Then  the  covariance  matrix  Pk  is  updated  as 

Pk={I~KkHk)Pk,k_ , 


(82) 


(83) 


(84) 
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The  process  then  repeats  starting  at  Equation  (80)  allowing  for  recursive,  on-line 
estimation  of  the  parameters  contained  in  x; . 

B.  PARAMETER  ESTIMATION  OF  A  DYNAMIC  SYSTEM 

The  system  identification  techniques  presented  to  date  are  developed  for  a  static 
system  of  form  yj  =  cp)d'  +(p]91  +...+(p"'9m  =$6 .  When  identifying  parameters  for  a 
dynamic  system  some  additional  steps  are  required  to  write  the  system  in  the  form  of  (46) 
[11].  Consider  the  dynamic  system: 

x  =  f0(x,u)  +  f(x,u) 

(ojj 

y  =  Cx 

0  is  a  matrix  of  unknown  parameters  and  O (x,u)  is  a  known  regressor  that  consists  of 
known  basis  functions.  This  dynamic  model  must  be  converted  to  a  static  system  in  order 
to  apply  the  system  identification  techniques  of  the  previous  section.  This  can  be 
accomplished  through  the  introduction  of  filtered  signals.  First,  rewrite  Equation  (85)  as: 

x  +  ax  =  ax  +  f0(x,u)  +  4>(x,  m)0  (86) 

A  filtered  version  of  x,  Xf,  is  also  introduced  as 

kf  +  axj-  =  ax  (87) 

where  a  is  the  filtering  time  constant.  Next,  define  z  =  x-xf ,  then  z  =  x-xf,  allowing  eq. 
(86)  to  be  rewritten  as: 

z  +  az  =  fQ(x,u)  +  <S>{x,u)®  (88) 

The  right  hand  terms  can  be  thought  of  as  forcing  functions,  and  the  first  order  ordinary 
differential  equation  has  a  known,  unique  solution: 

z  =  e‘"z(0)  +  je  T\f0(x,u)dr  +  \e  T^$>(x,u)dz®  (89) 

o  o 

Assume  the  initial  conditions  x(0)=x/0)  so  that  z(0)=0.  Then  equation  (88)  can  be 
simplified  to 

z(7)  =  ®0  +0^0  (90) 
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where  ®0  =je  a{t  T^f0(x,u)dr  is  the  filtered  version  of fo(x,u)  and  ®  ,  =  Je  a(t  r)<t>(x,u)dr 
0  0 

is  the  filtered  version  of  <P/x,u).  Furthermore,  rewriting  the  measurement  equation  in 

terms  of  z  and  x/,  the  dynamic  system  is  converted  into  a  static  system: 

v  -  Cxf  -  C®0  =  C®/©  (91) 


Since  &o  and  0/are  filtered  states,  they  can  be  calculated  recursively  as: 

+MX’U) 

®  f  =  -a®  f  +  f  (x,  u ) 


(92) 


Now,  the  dynamic  system  is  presented  as  an  equivalently  static  system  consisting 
of  filtered  version  of  the  known  and  unknown  dynamics  that  make  up  the  system. 
Because  of  this,  the  previously  derived  system  identification  techniques  can  be  applied  to 
the  equivalent  static  system. 


C.  PERSISTENCE  OF  EXCITATION 

The  quality,  or  richness,  of  the  signal  used  to  excite  a  system  during  parameter 
estimation  is  of  great  importance  for  the  quality  of  the  estimate.  For  the  GE  it  is  required 
to  ensure  convergence  to  the  true  values  of  the  parameters  being  estimated.  Recall  the 
equation  for  parameter  estimation  error  of  a  GE,  Equation  (66).  Solving  this  differential 
equation  gives 

i 

0,  =  0O  exp(-J  p0(pT  (r)tp(r)dr)  (93) 

0 

where  r  is  a  dummy  variable  used  for  integration  across  the  input  signal  (p(r) .  This 
implies  that 

i 

lim  [  cpT  ( r)<p{r)dr  =  go  (94) 

i— >oo  J 
0 

Therefore,  ©;.  will  converge  to  zero  (and  therefore  the  estimated  parameters  will 
converge  to  the  true  values)  if  there  are  positive  constants  T  and  a  such  that  for  all  /  >  0 
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(95) 


i+T 

J  <pT  (r)(p(r)dr  >  al 

o 


where  /  is  the  identity  matrix. 

Recall  from  Equations  (51)  that  to  guarantee  a  unique  solution  to  the  RLLS 
problem  <nro  must  not  be  singular.  This  concept  is  similarly  discussed  in  depth  by 
Astrom  in  [8]  as  requiring  the  matrix: 


2>2o--d  Z  u(i  -  Y)u(k  -  2)  -1  )u(k  —  n) 

72+1  22+1  72+1 

k  k  k 

2>(/-lM*-2)  Z»2(t-2)  ...£  u(i  -2)u(k  -n) 


72+1 


72+1 


72+1 


A. 

»(»-!)«(*- n) 


fu\k-n) 


72+1 


(96) 


to  be  of  full  rank  where  2/  denotes  the  input  signal,  i  is  the  time  step,  k  is  the  total  number 
of  time  steps,  and  n  denotes  the  order  of  the  system.  Astrom  defines  this  as  the  excitation 
condition.  For  long  data  sets  all  sums  in  Equation  (96)  are  taken  from  1  to  k  giving: 


C  =  lim  — ®r®  = 

k^co  k 


c(0) 

c(l) 


c(l) 

c(0) 


c(n  - 1)  c(n  -  2) 


c(n  - 1) 
c(n  -  2) 


c(0) 


(97) 


where  c(i)  are  the  empirically  detennined  covariances  of  the  input  signal  such  that 

1  k 

c(i)  =  \im  —  'Y'u(i)u(i-k)  (98) 

k 


Interpreting  Equations  (96)-(98)  helps  define  a  signal  as  PE  of  order  n  if  the  limit 
of  (97)  exists  and  Cn  is  positive  definite.  It  is  important  to  recall  here  that  the  only 
guarantee  required  for  existence  of  a  unique  solution  to  the  RLLS  problem  is  to  satisfy 
that  Equation  (51)  is  full  rank  which  is  primarily  driven  by  the  regressor  structure  in  ®  . 
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Therefore,  if  a  non-PE  signal  is  applied  to  a  system  the  RLLS  estimated  values  will 
converge  to  the  true  values,  but  only  slowly  Therefore,  if  it  is  possible,  a  PE  signal 
should  always  be  used  to  guarantee  rapid  convergence. 
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V.  SYSTEM  IDENTIFICATION  APPLIED  TO  VLBV300 


The  generic  hydrodynamic  model  for  an  underwater  vehicle  was  specialized  for  the 
vLBV300  by  introducing  various  assumptions  in  Section  III.B.  Even  with  these 
simplifying  assumptions  (often  times  ignoring  off-diagonal  coefficients),  the  resulting 
model  is  highly  coupled  (see  Equation  (46)  and  15  distinct  parameters  to  be  estimated). 
By  exciting  specific,  uncoupled  modes  for  the  vLBV300  sequentially,  individual 
parameters  can  be  isolated  and  thus  estimated.  For  example,  by  commanding  thrust  in 
the  surge  direction  through  Xprop,  only  motion  in  the  surge  direction  is  induced,  reducing 
Equation  (46)  to: 


Xu\u\u  11  +  X 


prop 


m  -  Xu 


v  =  0,  w  =  0,  p  =  0,  q  =  0,  r  =  0 


(99) 


The  filtering  techniques  of  the  previous  chapter  can  then  be  applied  to  this  simplified 
system  to  obtain  a  static  system  and  to  perform  the  system  identification: 


v,  =  % 


0  = 


m  -  X, 


1 


m  -  X, 


u  ill 


prop 


(100) 


The  inertia  parameters  are  not  identified  specifically  since  they  do  not  appear  in 
the  dynamic  equations  without  their  corresponding  added  mass  tenn  (e.g.,  Ix  and  K 

always  appear  together).  This  does  not  affect  the  applicability  of  the  mode  and  the 
assumed  values  presented  in  Chapter  III  are  used  for  clarification. 

Since  the  SeaBotix  ROV  is  also  controllable  in  sway,  heave,  roll  and  yaw,  the 
following  simplifications  are  similarly  possible  and  the  system  identification  of  the 
previous  section  can  similarly  be  applied. 

Sway: 
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h,„„W  +7, 


prop 


m  -  Y- 


it  =  0,  w  =  0,  p  =  0,  q  =  0,  r  =  0 


(101) 


Heave: 


Roll: 


Yaw: 


.  Z^w\WW+Zprop 

w  — - 

m  -  Zw 

it  =  0,  v  =  0,  p  =  0,  q  =  0,  r  =  0 


K 


P  = 


p\p\P  P  +  FBK  +  K prop 


[x-Kp 


it  =  0,  v  =  0,  w  =  0,  q  =  0,  r  =  0 


Nr\r\r  r|  +  ^ 


prop 


h~Nr 

it  =  0,  v  =  0,  w  =  0,  p  =  0,  q  =  0 


(102) 


(103) 


(104) 


However,  control  in  the  pitch  degree  of  freedom  is  not  possible  (Mprop= 0),  but  this 
mode  can  be  excited  by  taking  advantage  of  the  coupling  and  commanding  a  combination 
of  other  modes.  An  examination  of  the  dynamic  equation  associated  with  pitch  reveals 
that  by  exciting  the  vehicle  in  the  surge  and  heave  degrees  of  freedom,  the  pitch  portion 
of  Equation  (3 1)  can  be  excited,  resulting  in  the  simplified  to  v  =  o,p  =  0,r  =  0  and: 


u  = 


~(m  ~  Z-V)wq  +  Xmuu  +  X 
m~xu 

(m-X-)uq  +  ZMMw\w\  +  Z 


prop 


W  =  - 


prop 


q  =  - 


m~Zw 

~(ZW  -  Xu  )«W  +  Ma\qp  kl  +  FBM 


(105) 


(7y  —  M  ■  ) 
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In  terms  of  Equation  (85)  define: 


/n  = 


-{m  -  Z-v)wq  +  Xmu  \u\  +  X prop 

m~Xu 

0 m  -  X ■  )uq  +  Zwmw\w\  +  Zprop 

7  . 

‘‘w 

0 


m~zw 


(106) 


o  = 


0  0 

0  0 

?l?l  fbm-(Zw~Xu>w. 


0  = 


Mm  1 

(Iy-m.)  (. iy-m .) 


(107) 


(108) 


In  an  effort  to  verify  the  system  identification  approach  and  investigate 
persistence  of  excitation  for  the  system,  a  simulator  is  developed  according  to  the 
dynamics  presented  in  Equation  (46).  Chen  et  al.  [5]  presents  a  model  for  the  Seamor 
hovering-class  ROV  (the  parameters  were  estimated  using  the  Projective  Mapping 
Method,  while  the  mass  and  inertia  parameters  are  listed  as:  mass  =  20.4  kg,  Ix  =  0.429 
kgm  ,  and  Iy=  Iz=  0.609  kgm“).  The  results  described  in  [5]  are  used  as  true  values  and 
the  four  system  identification  techniques  described  in  Chapter  IV  are  in  turn  applied  to 
the  simulated  system  to  verity  that  the  true  parameters  are  estimated. 


A.  RLLS  ESTIMATION 

1.  RLLS  Estimation  Applied  to  Simulator 

The  hydrodynamic  coefficients  are  time-invariant  in  the  simulator  and  a  certain 
amount  of  measurement  and  process  noise  is  expected.  Recall  from  Equation  (51)  that 

<D  ®  must  be  of  full  rank  to  guarantee  the  RLLS  filter  will  converge  to  the  true  values  of 
the  parameters  being  estimated.  As  an  example,  the  surge  direction  Equation  (99)  is 
investigated  in  this  context.  Here  <p.  =  [«|«|  xprop^  so 

i= 1 


(«ihl)2 

u\u\X 

i  \  i  \  prop,i 

X  U  M- 

prop, r  i  \  i  | 

X2 

prop,i 

(109) 
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The  rank  of  the  matrix  described  by  Equation  (109)  was  analyzed  using  data  from  the 
previously  developed  simulator  and  determined  to  be  2,  or  full  rank.  Therefore, 
according  to  the  construction  of  the  RLLS  filter  for  this  system  the  estimated  parameter 
values  will  converge  to  the  true  values,  given  enough  data  (and  time).  This  result  is 
typical  for  the  other  DOF  equations  as  well. 

The  RLLS  estimator  is  applied  to  the  simulator  (with  coefficients  presented  in  [5]) 
to  sequentially  identify  the  hydrodynamic  coefficients  and  verify  the  correct 
implementation  of  the  RLLS  method.  As  an  example,  results  of  the  simulator  RLLS 
surge  testing  are  presented  in  Figure  13. 
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Figure  9.  Convergence  of  simulator  surge  parameters 


As  discussed  in  Section  IV. A  and  in  [8],  persistence  of  excitation  is  a  construct 
for  guaranteeing  speed  of  convergence  of  the  solution  to  the  RLLS  problem.  The 
estimated  parameter  values  can  be  seen  to  converge  to  the  true  values  rapidly,  implying 
that  the  system  is  persistently  excited  with  a  step  function.  A  “stair-step”  series  of  step 
inputs  are  more  suitable  if,  as  in  this  case,  testing  space  and  time  are  limited  (due  to  the 
test  tank  facility). 
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These  results  for  the  surge  direction  are  typical  for  the  other  coefficients  with  the 
exception  of  pitch.  This  is  due  to  un-modeled  dynamics  represented  by  FBm  in  Equation 
(37).  Chen  et  al.  [5]  did  not  present  the  value  for  zB,  the  position  of  the  center  of 
buoyancy  relative  to  the  center  of  gravity.  Without  knowing  this  value,  assumed  to  be 
zero  in  the  simulation,  exact  results  cannot  be  achieved.  The  restorative  force  due  to  zB 
shows  up  in  the  dynamics  for  the  roll  DOF  as  well  but  the  contribution  of  zB  to  FBk  is 
small  due  to  small  induced  pitch  angles  and  as  such  the  parameters  were  still  estimated 
accurately.  Full  estimation  results  are  presented  in  Table  3. 


True  value 

Est.  value 

Units 

X. 

-27.08 

-27.08 

kg 

Xu\u\ 

-61.117 

-61.117 

kg/m 

Y. 

-25.952 

-25.952 

kg 

kh 

-139.81 

-139.81 

kg/m 

z. 

-68.576 

-68.576 

kg 

z , 

w|  w| 

-51.724 

-51.724 

kg/m 

kp 

-61.683 

-61.683 

kgnr/rad 

K  ,  , 
p\p\ 

-12.0 

-12.0 

kgnrVrad2 

M * 

-79.411 

-82.363 

kgnrVrad 

M  , , 

?|?| 

-56.61 

-58.77 

kgnrVrad2 

N, 

-0.154 

-0.154 

kgnrVrad 

-1.772 

-1.772 

kgnrVrad2 

Table  3.  Comparison  of  hydrodynamic  parameters 


From  these  results,  it  can  be  concluded  that  the  RLLS  estimator  performs  properly  and 
that  the  signal  used  to  excite  the  system  contains  sufficient  PE. 

2.  RLLS  Applied  to  SeaBotix  vLBV300 

Since  the  vLBV300  is  not  yet  equipped  with  an  onboard  Inertial  Navigation 


System  (INS),  an  external  motion  capture  system  (VICON)  is  used  to  measure  the  vehicle 
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position  and  orientation  in  the  NPS  CAVR  test  tank.  From  these  pose  measurements,  the 
linear  and  angular  velocities  are  estimated.  The  VICON  system  consists  of  Infrared  LED 
arrays  and  cameras  that  track  reflective  markers  in  the  operating  space.  Due  to  the 
absorption  of  electromagnetic  signals  in  water,  the  current  system  is  only  applicable  for 
in-air  operation.  To  overcome  this  limitation,  the  vLBV300  has  been  extended  above  the 
water  surface  with  a  light-weight,  low  inertia  structure.  By  tracking  this  structure  and 
performing  the  appropriate  coordinate  transformations,  the  submerged  vehicle’s  motion 
can  be  tracked.  This  structure  has  a  small  effect  on  the  dynamics  of  the  vehicle,  and  a 
more  appropriate  solution  is  being  investigated  (such  as  relying  on  INS  data  instead). 
The  VICON  system  provides  high-accuracy  data  (<lcm)  at  high  data  rates  (100  Hz). 


Figure  10.  SeaBotix  vLBV300  in  the  instrumented  NPS  dive  tank 


There  are  some  limitations  associated  with  this  experimental  setup,  the  primary  of 
which  is  the  size  of  the  NPS  dive  tank.  At  approximately  1.4  meters  deep,  4.5  meters 
wide  and  6.5  meters  long  the  size  of  the  tank  prohibits  extended  data  collection,  an  effect 
that  is  exacerbated  during  higher  speed  runs.  In  order  for  the  VICON  system  to  provide, 
data  it  must  be  able  to  see  the  motion  capture  “pucs”  (IR  reflectors)  on  top  of  the 
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vehicle’s  extension.  As  the  vehicle  approaches  the  extremes  of  the  tank,  these  reflectors 
may  become  occluded  from  the  cameras,  causing  jumps  in  the  data  and  further  limiting 
the  space  available  for  experiments.  Another  limitation  is  the  VICON  system  itself. 
While  VICON  provides  very  high  accuracy  in  position  and  pose  data  it  is  a  physical 
system  and  is  therefore  subjected  to  some  measurement  noise.  VICON  does  not  provide 
velocities  or  angular  rates  and  therefore  the  provided  position  data  must  be  filtered  and 
differentiated  (and  transformed  as  required)  to  obtain  the  required  body  velocities  and 
angular  rates.  This  differentiation  has  the  effect  of  magnifying  the  noise  present  in  the 
VICON  measurements.  A  low-pass  filter  is  used  to  remove  high  frequency  noise  in  the 
measurements  but  those  effects  are  still  apparent  in  the  output  velocities  and  angular 
rates.  In  the  near  future,  CAVR  will  be  installing  an  INS  into  the  SeaBotix  vLBV300. 
This  will  help  to  improve  the  parameter  estimation  results  by  removing  the  need  to 
externally  measure  and  process  position  data  as  well  as  allowing  sustained,  at-sea  data 
runs.  All  of  these  effects  are  seen  in  the  application  of  RLLS  to  the  SeaBotix  vLBV300 
but  do  not  prevent  system  identification  results  from  being  obtained. 

One  of  the  fundamental  assumptions  for  the  generic  hydrodynamic  coefficients 
for  the  system  introduced  in  Chapter  III  is  time  and  state  independence.  This  may  be  a 
valid  assumption  for  systems  operating  around  a  trim  point,  as  is  implicitly  the  case  for 
the  simulator  (and  thus  the  results  for  the  previous  section).  However,  in  practice  these 
coefficients  are  dependent  on  the  system  state,  in  particular  vehicle  velocity.  This  can  be 
seen  in  the  initial  application  of  the  RLLS  estimator  to  the  vLBV300  using  real-world 
data.  The  RLLS  estimator  is  used  to  estimate  the  surge  parameters  of  the  vLBV300 
during  a  run  with  a  PWM  command  of  20  to  the  aft  thrusters  only.  The  resulting 
coefficients  are  used  to  obtain  a  response  which  is  compared  to  a  second  experimental 
run  also  conducted  at  a  PWM  of  20  to  the  aft  thrusters.  The  surge  displacement  results 
are  presented  in  Figure  1 1  and  the  velocity  results  are  presented  in  Figure  12. 

The  simulated  results  did  not  match  the  experimental  results  as  expected,  and  in 
fact  differed  by  a  significant  amount.  The  discrepancy  derives  from  the  speeds  at  which 
the  coefficients  were  derived  and  the  speed  at  which  the  vehicle  moved  in  the 
experimental  data.  The  average  steady  state  surge  velocity  u  of  the  run  used  to  determine 
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the  coefficients  was  0.3635  m/s  while  the  value  of  u  for  the  experimental  run  used  for 
comparison  was  0.2915  m/s.  The  effects  of  this  are  clearly  seen  in  Figures  1 1  and  12. 


Surge  Displacement  Comparison 


Figure  11.  Simulated  surge  displacement  compared  to  experimental  surge  displacement 

at  a  PWM  command  of  20 


Surge  Velocity  Comparison 


Figure  12.  Simulated  surge  velocity  compared  to  experimental  surge  velocity  at  a  PWM 

command  of  20 
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One  approach  to  capture  the  dependence  of  these  coefficients  on  velocity  in 
particular  is  to  introduce  dimensionless  coefficients.  Fossen  [1]  presents  a  table  of 
normalization  variables  (the  Prime  I  system),  which  is  used  to  achieve  the  desired  non- 
dimensional  form.  The  nonnalization  variables  relevant  to  this  work  are  presented  in 


Table  4. 


Units 

Prime-system  I 

Mass 

2 

Inertia 

2 

Force 

^ U2L2 

2 

Moment 

?U2L3 

2 

Quadratic  Damping 

?L2 

2 

Linear  Damping 

—UL2 

2 

Table  4.  Normalization  variables  used  in  the  Prime  I  system  (From  [1]) 


L  is  the  characteristic  length  of  the  vehicle  which  defines  the  scale  of  the  vehicle. 
In  this  case  it  is  chosen  to  be  the  length  of  the  vehicle  along  its  x  axis,  or  0.625  meters. 
p  is  the  density  of  water  at  the  temperature  of  the  tank  in  which  the  vehicle  is  being 
tested.  For  this  work  that  is  assumed  to  be  the  density  of  freshwater  at  15  degrees 
kg 

centigrade,  p  =  1,000 — - .  Finally,  let  the  steady  state  velocity  of  the  vehicle  be 
m 

U  =  y/u~  +  V"  +W~  (110) 

In  the  following,  the  subscript  c  applied  to  U  denotes  the  value  used  to  derive  the 
coefficients  while  the  subscript  r  applied  to  U  denotes  the  value  used  during  the 
comparison  run. 

Using  the  variables  from  Table  4,  it  is  possible  to  define  a  prime  set  of 
dimensionless  variables.  As  an  example,  the  variables  from  Equation  (99)  for  surge  are 
shown  in  detail  here  while  only  the  results  of  the  other  two  planar  DOF  are  presented. 
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Ill 


111  = 


(111) 


X,'  = 


x  ,;  =  ■ 


Pi? 

2 

X 


(112) 


X 


prop 


pl} 


V 

prop 

U2L2 


(113) 


(114) 


Inserting  Equations  (11 1)-(  1 14)  into  Equation  (99)  gives  the  non-dimensional  equation 

X,  '  tt2 


U  =  ■ 


"K|«l  I  I  . 

-u\u\  +  - 


L(m'-  Xu ')  1  1  Urn X u ') 


-X. 


(115) 


Similarly,  then 


Yw'  ,  ,  ry 

•  v  v  O'  , 

v  = - u - v  |v|  + 


Um'-YZ)  1  1  L(m’-Y7) 


-7. 


I  V  l\  Pr°P 


(116) 


w  =  - 


J  W  W  I  I  i 

-ww + 


L(m'-ZJ) 


I  -7  Pr°P 


Hm'-ZZ) 


(117) 


Since  Equation  (104)  involves  an  inertia  term,  I~,  instead  of  a  mass  term  as  before, 
the  coefficients  become 


v= 


PI? 


Nt '  = 


N, 

PL> 


(118) 


(119) 
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N'  =  - 


N, 


—1} 


N. 


N 


prop 


P  U2L3 


(120) 


(121) 


Inserting  Equations  (11 8)-(  121)  into  Equation  (104)  gives  the  non-dimensional  equation 


N, 


r  - 


u: 


L\IZ'~N 7)  1  1  L  (I z  '-N-') 


-N 


prop 


(122) 


It  is  easy  to  see  that  Equations  (11 5)-(  117)  and  (122)  can  still  be  re-arranged  to  fit 
the  form  of  Equation  (47)  so  the  RLLS  may  still  be  used  to  estimate  the  unknown 
parameters.  To  fully  investigate  the  vehicle  and  verify  the  normalization  scheme  used, 
two  surge  runs  are  presented  at  high  speed  (PWM  command  of  50)  and  low  speed  (PWM 
command  of  20)  to  the  aft  thruster  pair.  Dimensionalized  and  non-dimensionalized 
coefficients  are  determined  for  both  runs  and  initially  compared  to  different  runs  at  the 
same  PWM  command  that  the  coefficients  are  estimated  at.  Slightly  different  steady 
state  velocities  are  determined  for  the  two  runs  at  the  same  PWM  command  for  both  the 
high  and  low  speed  runs.  The  high  speed  run  results  are  presented  in  Figure  13  and  the 
slow  speed  results  are  presented  in  Figure  14. 
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Surge  Displacement  Comparison 


Figure  13.  Simulated  surge  displacement  results  compared  to  measured  surge 

displacement  results  at  a  PWM  command  of  50  using  coefficients  determined 

at  a  PWM  command  of  50 


Surge  Displacement  Comparison 


Figure  14.  Simulated  surge  displacement  results  compared  to  measured  surge 

displacement  results  at  a  PWM  command  of  20  using  coefficients  determined 

at  a  PWM  command  of  20 


54 


In  both  cases  it  can  be  seen  that  the  non-dimensional  simulator  performs  more 
accurately  than  the  dimensionalized  version,  but  in  the  high  speed  run  the 
dimensionalized  simulator  also  performs  well.  Additionally,  when  coefficients  derived 
during  the  low  speed  run  are  used  to  simulate  a  high  speed  run  (or  vice  versa)  the  results 
diverge  significantly  from  the  experimental  data.  This  is  in  part  due  to  the  application  of 
the  model  away  from  the  trim  condition.  Figure  15  shows  the  divergence  of  the  high 
speed  coefficients  used  to  simulate  a  low  speed  run.  Low  speed  coefficients  simulating  a 
high  speed  result  also  demonstrate  this  effect  and  diverge  significantly. 


Surge  Displacement  Comparison 


Figure  15.  Comparison  of  surge  displacement  to  experimental  displacement  using 
coefficients  determined  at  PWM  command  of  50  and  a  simulated  run  at  a 

PWM  command  of  20. 


Both  the  inaccuracy  in  the  dimensionalized  simulator  seen  in  Figures  13  and  14 
and  the  divergence  of  the  non-dimensionalized  simulator  shown  in  Figure  15  can  be 
explained  by  an  examination  of  the  Reynolds  number.  The  Reynolds  number  is  a 
dimensionless  parameter  that  relates  the  viscous  behavior  of  all  Newtonian  fluids  [20].  It 
is  expressed  mathematically  as 


Re  = 


pVL 


A 


(123) 
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where  V  and  L  are  the  characteristic  velocity  and  length  that  describe  the  vehicle’s 
passage  through  the  fluid,  p  is  the  density  of  the  fluid  the  vehicle  is  travelling  in,  and  p 
is  the  viscosity  of  the  fluid.  The  Reynolds  number  is  also  used  to  determine  if  a  condition 
of  similarity  exists  between  the  results  detennined  from  two  different  experiments. 
White  defines  the  similarity  condition  as  existing  if  the  Reynolds  number  of  one  run  is 
equal  the  Reynolds  number  of  the  second  run.  He  demonstrates  this  using  a  ratio  of 
generic  force  coefficients  [20] 


A 

f2 


A 
P  2 


V2 

\V2J 


V 

l2 


(124) 


According  to  Equation  (124)  a  scaled  relationship  between  the  two  force  coefficients 
exists  if  the  Reynolds  number  of  the  two  experiments  is  the  same.  For  this  work,  since 
the  density  of  the  water  and  the  characteristic  length  of  the  vehicle  is  the  same  between 
experiments  Equation  (124)  simplifies  to 


E 

E 


(125) 


which  shows  that  the  velocities  must  be  equal  in  order  for  the  condition  of  similarity  to 
exist.  Because  the  experiments  are  conducted  at  different  velocities  the  Reynolds 
numbers  are  not  the  same  and  the  coefficients  cannot  be  scaled.  This  difference  is  not 
noticeable  in  Figures  13  and  14  because  the  steady  state  velocities  are  similar  enough, 
(i.e.,  it  is  close  enough  to  the  trim  condition).  This  is  not  the  case  for  the  high-vs  low- 
speed  runs. 


The  second  conclusion  drawn  from  an  analysis  of  the  Reynolds  number  explains 
the  difference  between  the  dimensionalized  and  non-dimensionalized  simulators  in  the 
low  speed  analysis.  Recall  another  assumption  previously  made  that  only  turbulent  flow 
would  exist  around  the  vehicle  during  its  motion.  This  allowed  for  simplification  of  the 
damping  matrix  D(v )  to  only  include  the  quadratic  damping  terms  and  allowed  for  the 
neglect  of  the  linear  damping  terms. 
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The  Reynolds  number  is  also  a  good  indicator  of  which  type  of  flow  will  be 
present  around  a  body  moving  in  a  fluid.  White  states  that  turbulent  flow  will  be  present 
at  Reynolds  numbers  greater  than  approximately  106  [20],  Calculating  the  Reynolds 
numbers  for  the  low  speed  and  high  speed  runs  gives 

1000^0.3635  —  0.625kg 

Re.  = - ^ - =  2.06*1 05 

,ow  1. 1*1 0'3  Pa  -s 

1000^-0.8630  —  0.625kg 
Re...  = - m - - =  4.903*1 05 


"'high 


1.1*10  Pa-s 


(126) 


These  Reynolds  numbers  show  the  flow  around  the  vehicle  should  be  totally 
laminar  for  the  low  speed  run,  not  turbulent,  and  therefore  that  ignoring  the  linear 
damping  terms  is  not  a  good  assumption.  This  effect  is  particularly  visible  in  the  low 
speed  run  shown  in  Figure  14.  The  fact  that  the  high  speed  run  in  Figure  13  matches  well 
implies  that  for  the  vLBV300  turbulent  flow  is  beginning  to  dominate  around  0.9  meters 
per  second  therefore  the  vehicle  is  likely  operating  at  the  end  of  the  transition  region. 
The  damping  matrix  of  the  model,  then,  must  be  updated  to  include  the  linear  damping 
coefficients  instead  of  quadratic  damping  for  the  low  speed  run. 

Dip)  =  -diag  ([Xu ,  Yv ,  Zw,Kp  ,Mq,Nr])  (1 27) 


Equations  (99)-(  104)  then  become 


11  = 


Xuu  +  X  prop 

m  -  Xu 


and  Equation  (108)  becomes 


V  +  Yprop 
m  -  Y- 

ZwW  +  Zprop 
m  -  Z- 


(128) 


P  = 


KpP  +  FBK  +  K  prop 


'x-Kp 


N/  +  Nprop 

IZ-N, 
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-(m-zw)wq  +  xuu  +  x  pmp 
m  ~  Xu 

(m  -  Xu  )uq  +  zww  +  z prop 
m~Zw 


(129) 


q  = 


-(Zw-Xu>w+Mqq+FBM 


Hy-m.) 


RLLS  regression  is  again  applied  to  the  data  sets  surge,  sway,  heave,  and  yaw 
using  the  updated  model  regressor  form 


u,  = 


1 

u 

1 

1 

£ 

- 1 

i 

V 

prop 

(130) 


as  well  as  the  updated  non-dimensionalized  regressor  form.  The  surge  equations  are 
presented  as  an  example  because  the  other  updated  DOF  equations  are  typical  to  surge. 


yt  = 


1 

c: 

c: 

to 

1 _ 1 

u 

Lim'-XJ  L(m'-X-) 

v' 

prop 

(131) 


X  =■ 


pn 


(132) 


The  surge  displacement  results  using  the  linear  damping  matrix  for  a  slow  speed 
run  are  presented  in  Figure  16  and  show  that  the  damping  model  is  still  not  entirely 
correct. 
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4.5 


Surge  Displacement  Comparison 


Figure  16.  Surge  displacement  comparison  using  linear  only  damping  model 

Based  on  these  results  it  is  concluded  that  the  true  force  due  to  damping  is 
probably  a  combination  between  linear  and  quadratic,  between  laminar  and  turbulent. 
While  the  Reynolds  number  for  this  run  suggests  total  laminar  flow,  the  equation  used  for 
this  is  the  Reynolds  number  of  an  infinitely  flat,  smooth  plate.  The  SeaBotix  vLBV300  is 
clearly  not  a  smooth,  flat  plate.  Thus,  it  is  more  likely  that  the  vehicle  operates  in  the 
transition  region  from  laminar  to  turbulent  flow.  Most  likely,  there  are  local  instances  of 
turbulent  flow  due  to  the  rough  surface  of  the  sides  and  underside  of  the  vehicle  while 
laminar  flow  exits  on  the  largely  flat  upper  surface.  The  damping  matrix  can  be 
represented  as  a  combination  of  both  the  linear  and  quadratic  terms  and  this  is  a  good 
area  for  further  work. 

Since  the  quadratic  damping  model  performed  well  at  a  PWM  command  of  50  for 
both  the  dimensionalized  and  non-dimensionalized  simulators  the  full  results  for  the  high 
speed  experiments  in  surge,  sway,  heave,  and  yaw  are  presented  here  as  a  starting  point 
for  future  work.  As  demonstrated,  these  values  will  only  be  accurate  around  the  speed  at 
which  they  are  derived  but  they  are  useful  as  a  starting  point  for  future  work. 
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Complete  system  identification  results  for  the  SeaBotix  vLBV300  at  high  speeds 
are  presented  in  Table  5.  The  high  speed  runs  for  surge  and  sway  are  conducted  at  a 
PWM  command  of  50  while  the  heave  and  yaw  runs  were  conducted  at  PWM  command 
of  40  (the  shallowness  of  the  tank  in  the  heave  DOF  and  the  tether  wrapping  around  the 
vehicle  limited  operations).  The  convergence  of  the  surge  coefficient  values,  both 
dimensionalized  and  non-dimensionalized,  is  presented  in  Figure  17  and  the  sway  and 
yaw  DOF  results  are  similar  to  this.  The  heave  DOF  results  are  presented  in  Figure  18  to 
demonstrate  the  limitation  imposed  by  the  shallow  depth  of  the  tank.  In  order  to  prevent 
the  vehicle  from  colliding  with  the  bottom  of  the  tank,  only  about  two  seconds  of  data  is 
available  to  the  estimator.  Figure  18  shows  that  the  RLLS  appears  to  converge  then 
makes  a  correction  (both  in  the  dimensional  and  non-dimensional  filter)  at  the  very  end. 
Clearly,  longer  data  runs  would  provide  more  time  for  the  filter  to  converge  and  therefore 
more  accurate  results. 


Dimensional  Surge  Coefficient  Convergence 


Figure  17.  Convergence  of  surge  coefficients  at  high  speed 
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Dimensional  Heave  Coefficient  Convergence 


Non  Dimensional  Heave  Coefficient  Convergence 


Figure  18.  Convergence  of  heave  coefficients  at  high  speed 


Dimensionalized 

Units 

Non- 

Dimensionalized 

X. 

-13.5778 

kg 

-0.1146 

X  , 

u  u 

-27.7411 

kg/m 

-0.1416 

Y. 

-27.9347 

kg 

-0.2334 

Y  | 

V  v| 

-50.6868 

kg/m 

-0.2592 

zw 

-46.3258 

kg 

-0.3971 

-64.5970 

kg/m 

-0.3000 

Nr 

-3.1023 

kgnr/rad 

-0.8112 

N  | 

r  r 

-2.1709 

kgnr/rad2 

-0.3416 

Table  5.  Parameter  estimation  results  for  SeaBotix  vLBV300  at  high  speed 


3.  Conclusion:  RLLS  Estimator 

The  RLLS  estimator  is  applied  to  both  a  dimensionalized  and  non- 
dimensionalized  model  for  the  SeaBotix  vLBV300  platform.  The  benefits  of  the  non- 
dimensionalized  approach  are  demonstrated  to  make  small  adjustments  around  the  trim 
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conditions.  The  non-dimensional  RLLS  produced  both  a  low-speed  and  a  high-speed 
model  that  performed  well  when  tested  at  their  respective  speeds  and  generally  performed 
better  than  the  dimensionalized  simulator.  A  Reynolds  number  analysis  and  supporting 
experimental  data  suggests  that  the  vehicle  operates  in  different  regimes  for  the  low-and 
high-speed  operations  and  an  alternative  model  for  low-speed  operations  is  investigated. 
From  these  results  it  is  concluded  that  a  combination  of  these  models  may  be  required. 
Also,  since  there  is  a  strong  dependence  of  the  model  coefficients  on  vehicle  state 
(velocity  in  particular),  it  may  be  more  appropriate  to  model  the  system  with  time- 
varying  parameters,  which  can  also  be  estimated  with  system  identification  techniques 
(recommended  future  work). 

The  RLLS  is  an  excellent  estimator  for  this  work.  It  has  the  ability  to  estimate 
parameters  of  both  constant  and  time-varying  parameters,  which  as  demonstrated,  may  be 
a  useful  extension  of  this  work  to  provide  the  most  accurate  hydrodynamic  parameter 
results  possible.  A  major  current  limitation  is  the  size  of  the  NPS  dive  tank,  which  does 
not  allow  for  data  runs  longer  than  approximately  12  to  15  seconds  at  low  speed  or  five- 
to-seven  second  runs  at  high  speed.  In  the  presence  of  noise  and  uncertainty  this  is  often 
not  enough  time  for  the  filter  to  converge.  Longer  data  runs  would  greatly  benefit  the 
RLLS  by  allowing  it  more  time  to  converge. 

Additionally,  this  portion  of  the  work  identified  the  ideal  sequence  of 
commanded,  decoupled  motion  to  identify  the  full  set  of  model  parameters.  Because  of 
the  ability  to  perform  uncoupled  motion  in  the  surge,  sway,  heave,  roll,  and  yaw  DOF  the 
parameters  of  those  equations  should  be  determined  in  that  order.  This  is  required  to 
allow  for  the  identification  of  the  pitch  coefficients  through  coupling  of  surge  and  heave. 
After  these  primary  terms  have  been  detennined,  they  can  be  plugged  back  into  the 
appropriate  equations  of  motions  (used  in  subsequent  estimation  runs). 

This  strategy  can  also  be  applied  through  careful  analysis  of  the  equations  of 

motion  to  identify  more  coupling  terms  (to  relax  some  of  the  assumptions).  The  coupled 

motions  can  be  perfonned  in  all  possible  permutations  (surge-sway,  heave-roll,  sway- 

heave,  etc.)  in  order  to  determine  the  cross-coupling  terms  that  were  initially  neglected  in 

the  equations  of  motions.  Finally,  with  these  initial  parameters  in  place,  the  filter  can  be 
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deployed  on-line  on  the  vehicle  to  constantly  learn  and  update  from  these  baseline  values. 
As  more  data  is  collected  through  normal  operations  the  estimates  can  be  refined  or  can 
adapt  to  mid-mission  configuration  changes. 

B.  NEURAL  NETWORK  ESTIMATION 

1.  Neural  Network  Applied  to  Simulator 

The  strength  of  the  neural  network  approach  is  that,  using  sufficient  number  of 
nodes,  an  arbitrary  number  of  inputs  can  be  mapped  to  an  arbitrary  number  of  outputs 
with  good  accuracy  [10].  The  type  of  NN  investigated  is  the  Non-linear  Autoregressive 
(with  External  Input)  NN  (NARX).  This  type  of  model  architecture  is  generally  more 
accurate  than  other  models  because  it  incorporates  feedback.  That  is,  previous  values  of 
the  output  y  are  fed  back  to  help  improve  the  accuracy  of  the  network  and  continually 
train  the  network.  As  an  example,  thruster  force  in  the  surge  direction  was  mapped  to 
surge  velocity  for  the  simulated  Seamor  ROV.  Both  the  Xprop  command  and  the  resulting 
surge  velocity  are  corrupted  by  low  power,  white  noise.  First,  a  network  using  three 
nodes  and  two  delays  is  trained  then  a  network  using  10  nodes  and  two  delays  is  trained 
The  diagnostic  graph  of  the  NARX  response  for  the  three  node  network  is  presented  in 
Figure  19  while  the  response  for  the  10  node  network  is  presented  in  Figure  20. 

In  the  learning  block,  70  percent  of  the  supplied  data  was  used  to  create  and  train 
the  network,  15  percent  was  used  to  validate  the  network,  and  the  final  15  percent  was 
used  to  further  test  the  network.  The  top  graph  shows  the  response  of  the  output  element 
one,  or  for  this  network,  surge  velocity  u,  over  time  plotted  against  test  and  validation 
points.  The  lower  graph  shows  the  magnitude  of  the  error  between  the  network  output 
and  the  supplied,  true  output  data.  . 
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Response  of  Output  Element  1  for  Time-Series  1 
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Figure  19.  Diagnostic  plot  of  NARX  with  three  nodes  and  two  time  delays 


Response  of  Output  Element  1  for  Time-Series  1 
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Figure  20.  Diagnostic  plot  of  NARX  with  10  node  and  two  time  delays 
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A  comparison  between  Figures  19  and  20  shows  that  increasing  the  number  of 
nodes  from  three  to  10  has  little  effect  on  the  magnitude  of  the  error.  However,  a 
comparison  of  their  outputs  when  plotted  against  the  simulated  surge  velocity 
demonstrates  another  constraint  of  NN:  using  too  many  nodes  to  model  a  simple  system. 
This  comparison  of  the  surge  velocities  produced  by  both  networks  is  presented  in  Figure 
21. 


Comparison  of  Neural  Network  and  Measured  Response 


Figure  21.  Comparison  of  actual  surge  velocities  to  10  node  and  three  node  NARX 

Figure  2 1  shows  that,  for  the  simple  input-output  relationship  between  Xprop  and 
surge  velocity  a  three  node  network  is  capable  of  modeling  the  response  of  this  system. 
It  also  does  a  good  job  of  rejecting  the  high  frequency  noise  used  to  corrupt  the  inputs 
and  outputs.  The  10  node  network,  however,  attempts  to  model  the  unwanted  high 
frequency  noise  content.  This  demonstrates  the  power  of  increasing  the  size  of  a  NN  as 
well  as  the  dangers  of  doing  so.  Therefore,  when  choosing  a  NN  size  careful 
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consideration  must  be  given  to  using  enough  nodes  to  adequately  model  the  system  while 
ensuring  the  noise  inherent  to  the  system  is  not  also  modeled  by  the  NN. 

Another  important  consideration  when  modeling  using  NN  is  the  data  used  to 
train  the  network.  The  network  will  only  be  as  good  as  the  information  contained  in  the 
training  data.  For  example,  when  initially  training  the  NN  used  to  map  the  Xprop  to  surge 
velocity,  a  disproportionate  amount  of  the  data  used  was  for  the  steady  state  condition  of 
the  simulator  after  it  reached  cruising  speed.  Since  the  NN  toolbox  randomly  chooses 
points  from  the  data  provided  to  train  the  network,  this  resulted  in  poor  modeling  of  the 
transient.  Truncating  the  data  used  to  train  the  NN  to  involve  an  equal  amount  of  data  for 
the  transient  and  steady  state  generated  the  much  better  results  seen  in  Figure  21. 
Therefore,  if  a  single  NN  is  expected  to  model  an  entire  complex  system  then  when 
training  it  is  crucial  to  do  so  values  for  all  possible  operational  ranges  and  conditions  and 
to  ensure  all  dynamics  are  equally  represented  in  the  amount  of  data  provided  to  train  the 
system. 

Next,  a  10  node,  four  time  delay  NARX  neural  network  is  created  and  trained  for 
all  four  controllable  DOF  of  the  Seamor  ROV.  This  is  accomplished  by  exciting  the  four 
controllable  DOF  sequentially  in  the  order  surge,  sway,  heave,  yaw  then  exciting  the 
vehicle  in  the  coupled  modes  of  surge-sway,  surge-heave,  surge-yaw,  sway-heave,  sway- 
yaw,  heave-yaw,  surge-sway-heave,  and  surge-sway-yaw.  The  goal  of  this  input 
sequence  is  to  excite  the  vehicle  as  completely  in  all  six  DOF  (relying  on  coupling  effects 
for  roll  and  pitch)  as  possible  in  order  to  ensure  the  information  used  to  train  the  NARX 
is  a  rich  as  possible.  The  output  values  used  to  train  the  NARX  are  the  simulated  body 
velocity  and  angular  rate  values  u,  v,  w,  and  r.  The  input  commands  and  output 
measurements  are  corrupted  with  low  power,  white  noise  to  make  the  simulation  as 
realistic  as  possible.  The  diagnostic  results  relevant  to  the  surge  output  element  are 
presented  in  Figure  22. 

Because  of  the  computational  infrastructure  available  for  this  work,  trying  to 

include  roll  and  pitch  is  not  possible.  Nor  are  longer  data  runs.  The  computer  used 

cannot  handle  the  training  process  for  even  a  small,  10  node  network  using  six  inputs,  six 

outputs,  and  the  large  array  of  training  data.  It  is  also  difficult  to  increase  the  number  of 
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nodes  for  non-simple  (one  input  to  one  output)  systems  which  again  demonstrates  the 
importance  of  choosing  a  network  size  to  adequately  model  the  system  while  avoiding 
“over-modeling”  and  over-burdening  of  the  computer  used  to  create  and  train  the 
network. 


Response  of  Output  Element  1  for  Time-Series  1 
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Figure  22.  Diagnostic  results  for  surge  element  of  complex  NARX 


To  verify  this  model,  and  that  the  data  used  to  train  it  was  sufficient,  the  NARX 
was  tested  using  simple  coupled  motion  between  surge  and  sway  since  the  data  used  to 
train  this  NN  included  similar  coupled  motion.  The  results  of  plugging  the  resulting 
NARX  model  into  the  simulator  are  presented  in  Figure  23. 
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Figure  23.  Comparison  of  surge  velocity  response 


The  results  presented  in  Figure  23  show  the  power  of  using  information  rich 
training  data  and  choosing  the  right  network  size.  Because  a  complex  and  information 
rich  training  scheme  is  used  to  train  the  network,  when  a  different  input  signal  is  applied 
the  NARX  is  able  to  adapt  and  correctly  model  the  system  (with  some  difficulty  in  the 
transient  section  due  to  the  network  response  time).  This  serves  to  further  highlight  the 
importance  of  the  data  used  to  train  the  NN.  Also,  since  a  more  complex  relationship  is 
modeled,  a  greater  number  of  nodes  are  required  to  adequately  model  the  system.  . 

2.  Neural  Network  applied  to  SeaBotix  vLBV300 

Now  that  the  NN  concept,  specifically  using  the  NARX  architecture,  has  been 
validated  using  the  simulator,  it  is  applied  to  the  SeaBotix  vLBV300.  As  discussed  in 
Section  V.B,  while  limited  computing  power  restricts  the  number  of  nodes  that  can  be 
employed,  the  best  results  are  produced  ensuring  the  richest  possible  data  is  used  to  train 
the  networks..  The  NARX  experiment  applied  to  the  vLBV300  was  a  stair-step  input  the 
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surge  direction  where  only  the  thrust  generated  in  the  surge  direction  was  supplied  as 
input  and  only  u,  the  surge  velocity  was  used  as  the  targeted  output.  A  three  node 
network,  using  two  time  delays  was  created  in  MATLAB  with  the  error  analysis 
presented  in  Figure  24.  A  single  node,  single  delay  network  was  also  created  but  the 
diagnostic  results  are  not  presented  as  they  are  very  similar  to  the  three  node  network. 
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Response  of  Output  Element  1  for  Time-Series  1 
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Figure  24.  Diagnostic  results  for  velocity  mapping,  three  node,  two  delay  NARX  in  surge 

direction  only 


Again,  the  error  plot  demonstrates  just  how  well  a  very  small  NN  can  be  trained  to  follow 
a  non-linear  dynamic  system.  As  before,  the  network  was  then  verified  by  comparing  its 
outputs  to  measured  data  as  presented  in  Figure  25. 
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Comparison  ofNN  and  Experimental  Surge  Velocity 


Figure  25.  Comparison  of  three  and  one  node  NARXs  to  measured  data 

Even  the  single  node  NARX  was  able  to  match  the  transient — almost  too  well:  the  single 
node  network  is  able  to  match  the  noise  inherent  to  the  Vicon  system. 

Next,  a  complex  run  by  the  SeaBotix  coupling  the  surge,  sway,  heave,  and  yaw 
degrees  of  freedom  was  used  to  train  a  complex  NN.  Because  of  the  small  size  of  the 
tank  and  a  limitation  with  capturing  joystick  to  PWM  commands,  the  network  could  not 
be  sufficiently  trained  to  model  the  system.  Based  on  the  simulation  results,  this  is  due  to 
not  enough  data  being  collected  for  the  various  DOF.  The  simulation  results  are  achieved 
by  using  100  simulated  seconds  of  data  with  a  sample  rate  of  10  Hz.  This  yielded  over 
10,000  data  points.  It  was  not  possible  to  collect  this  much  information,  or  as  rich  as 


70 


information  using  the  actual  SeaBotix  due  to  the  current  experimental  setup  and  reliance 
on  externally  collected  motion  data. 

Conceptually,  a  single  well-trained  NN  will  be  able  to  approximate  motion  in  any 
DOF  or  number  of  coupled,  controllable  DOFs,  but  as  shown  that  will  require  a  network 
that  is  trained,  and  re -trained,  with  rich  and  comprehensive  input  data.  This  ideal 
network  may  require  a  large  number  of  nodes  depending  in  the  complexity  of  the  system 
being  modeled,  but  as  shown  too  many  nodes  will  actually  degrade  the  performance  of 
the  NN  estimator.  This  training  data  must  contain  information  on  all  operating  modes, 
controllable  DOF,  spanning  the  full  range  of  motions  of  the  anticipated  operations  in 
order  to  ensure  success.  This  highly  accurate  model  will  be  useful  for  both  motion 
prediction  and  fault  detection  and  assessment.  For  example,  if  the  vehicle  suspects  it  is 
damaged  due  to  onboard  diagnostics  or  because  of  divergence  of  true  motion  from 
another  model  (perhaps  supplied  by  the  RLLS  method)  it  could  enter  a  sub-routine  where 
it  performs  motion  similar  to  those  used  to  create  the  simplified  NN.  Significant 
deviation  from  this  simplified  motion  would  not  only  imply  damage  or  obstruction  but 
would  also  simplify  the  fault  diagnosis  process  since  fewer  numbers  of  thrusters  are  used 
as  inputs  to  these  networks. 

C.  GRADIENT  ESTIMATOR 

The  GE  is  greatly  limited  by  the  presence  of  un-modeled  disturbances  and 
measurement  noise,  and  is  very  sensitive  to  the  values  used  to  tune  the  filter  -  particularly 
step  size,  po  as  defined  in  Section  IV. A. 3.  As  in  Chapter  IV,  since  both  of  these  are 
expected  to  be  present  in  this  work  the  GE  is  not  a  good  choice  for  online  parameter 
estimation  of  a  ROV  and  is  not  further  investigated.  This  is  confirmed  by  attempting  to 
use  the  GE  to  estimate  the  surge  parameters  of  the  Seamus  ROV  presented  in  [5].  Recall 
that  the  only  guarantee  of  convergence  for  a  GE  is  sufficient  PE,  therefore  a  chirp  signal 
was  used  to  excite  the  simulated  system  since  it  a  signal  of  the  highest  order  of  PE  and  is 
rich  in  frequency  content  as  well. 
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First,  the  system  is  excited  using  a  chirp  signal  (a  sinusoid  based  signal  of  high  PE 
that  increases  its  frequency  over  a  set  period  time)  with  no  simulated  noise  and  the 
estimator  step  size,  po,  is  varied  from  one  to  1000.  The  results  presented  in  Figure  26  are 
typical  for  all  p0  values. 


Performance  of  Gradient  Estimator  for  the  Surge  Direction 


Figure  26.  GE  results  in  the  surge  direction  with  no  noise  added 

The  results,  while  oscillating,  can  be  seen  to  be  drifting  down  towards  the  true 
values  of  61.117  kg/nr  and  27.08  kg  for  X uu  and  X-,  respectively.  X u  even  arrives  in 

the  general  vicinity  of  its  true  value  and  then  proceeds  to  oscillate  around  it.  X u  u  , 

however,  is  very  far  from  the  true  value  even  after  a  30  second  period  of  excitation.  In 
order  for  the  Xu^  value  to  arrive  at  its  true  value  the  experiment  would  have  to  run  for  a 

very  long  period  of  time  (on  the  order  of  thousands  of  seconds  in  this  example)  Because 
of  the  small  size  of  the  NPS  tank  long  experimental  runs  are  not  possible.  As  such  a 
longer  simulated  experiment  was  not  conducted  since  an  operational  equivalent  is  not 
currently  possible  and  because  the  RLLS  estimator  is  shown  to  function  better  in  the 
presence  of  both  noise  and  uncertainty.  Combining  this  long  time  to  convergence,  with 
the  oscillation  caused  by  the  choice  of  po,  reveals  the  GE  to  be  a  poor  choice  for  this 
work. 
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Adding  even  very  low  power  noise  (a  magnitude  of  0.0001)  to  the  regressor 
inputs  to  the  filter  illustrates  the  vulnerability  of  the  GE  to  noise.  The  results  presented  in 
Figure  27  further  show  why  the  GE  is  a  poor  choice  for  this  work.  Because  of  this,  the 
GE  method  was  not  applied  to  the  SeaBotix  platform. 


Performance  of  Gradient  Estimator  for  the  Surge  Direction 


Figure  27.  GE  results  in  the  surge  direction  with  added  noise. 


D.  BAYESIAN  FILTERING 

By  choosing  a  zero  mean,  normal  distribution  for  the  probability  functions  related 
to  the  coefficients  to  be  estimated  as  well  as  the  measurement  and  process  noise 
distributions  the  Bayesian  filter  becomes  an  extended  Kalman  filter.  In  order  to 
successfully  apply  the  Kalman  filter  equations  a  system  must  be  observable  as  derived  by 
Ogata  [18]. 


To  check  for  system  observability,  first  the  system  must  be  linearized  and  written 
in  the  canonical  state-space  form 


x  =  Ax  +  Bu 
y  =  Cx  +  Du 


(133) 
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The  surge  direction  is  considered  first  by  re-arranging  Equation  (99)  into  the  form  of 
Equation  (133)  to  give 

u 
ax 
a2 

ii 

x=  ax 
_a2 

2  *ju,au  ut\u\  Xprop. 

A=  0  0  0  (134) 

0  0  0 

C  =  [  1  0  0] 

y-u 


In  Equation  (134)  ax  = 


(m~Xu) 


and  a2  = 


(m-X,) 


.  It  is  important  to  note  that  in 


Equation  (134)  u  refers  to  the  surge  velocity  in  the  body  frame  of  reference  and  in 
Equation  (133)  it  refers  to  a  control  input.  Matrices  B  and  D  are  ignored  here  since  they 
do  not  enter  into  the  observability  analysis.  Then,  from  [18],  observability  of  the  system 
is  given  by 


(135) 


where  n  is  the  size  of  the  state  space,  or  three,  in  this  example.  The  observability  matrix 
from  Equation  (135)  is  then  calculated  as 

■  c  1  r  i  o  o 

CA  =  2  ^au  ui\u\  Xpropi  (136) 

CA\  [(2^aX  i)2  2yJu'ialjUj  \ui \  aXiXpropi_ 


The  detenninant  of  matrix  A  is  zero  which  means  that  this  observability  matrix  is  rank 

deficient.  Therefore,  the  system  is  not  observable.  Since  the  system  is  not  observable  the 
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extended  Kalman  filter  cannot  be  applied  to  this  work.  This  is  not  an  atypical  result. 
While  most  physical  systems  are  in  fact  observable,  the  way  they  are  expressed 
mathematically  may  not  be  [18].  Therefore,  the  Bayesian  approach  is  not  investigated 
further  for  this  work  and  it  is  not  applied  to  the  SeaBotix  vLBV300. 
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VI.  CONCLUSIONS 


A.  SUMMARY 

Estimating  the  hydrodynamic  coefficients  of  an  underwater  vehicle  is  a  difficult 
process.  The  motion  of  a  body  in  areal  fluid  is  a  highly  complex,  highly  coupled,  and 
non-linear  process.  A  propulsion  model  is  developed  that  maps  individual  thruster  low 
level  PWM  commands  to  generated  thrust  and  a  geometric  model  that  converts  thrust  into 
body  forces  and  moments.  Then,  starting  with  the  first  principles  for  the  equations  of 
motion  for  a  six  DOF  submerged  vehicle  as  presented  by  Fossen  [1],  several  assumptions 
are  introduced  to  simplify  this  model  and  create  a  parametric  representation  of  the 
SeaBotix  vFBV300  THAUS. 

This  thesis  is  focused  on  applying  various  system  identification  techniques  to 
learn  the  model  parameters  for  the  THAUS.  Since  these  system  identification  techniques 
are  developed  for  static  regression  models,  a  method  of  converting  a  dynamic  system  into 
an  equivalent  static  system  is  presented  before  four  system  identification  techniques  are 
applied:  recursive  linear  least  squares  (RTFS),  computational  neural  networks  (NN), 
gradient  estimator  (GE),  and  a  Bayesian  filtering  method.  The  GE  and  Bayesian 
estimator  approaches  are  not  suitable  for  this  work  or  not  applicable  and  as  such  no 
results  are  compared  to  the  other  parametric  model  based  system  identification  method  - 
the  REES. 

The  RLLS  approach  is  applied  to  learn  the  coefficients  at  a  low-  and  high-speed 
trim  condition,  respectively.  By  sequentially  exciting  the  system  in  individual,  decoupled 
directions,  the  model  parameters  can  be  estimated.  Since  the  vehicle  cannot  be  controlled 
in  the  pitch  direction,  coupled  motion  between  the  surge  and  heave  directions  are 
required  for  system  identification.  It  is  shown  that  these  parameters  can  indeed  be 
estimated  online.  This  is  demonstrated  on  a  simulator  with  known  parameters  as  well  as 
the  vLBV300  platform,  utilizing  an  external  motion  capture  system  and  performing  tests 
in  the  NPS  CAVR  test  tank.  Model  accuracy  is  good  in  most  directions;  however  the 
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experimental  setup  limited  the  lengths  of  the  data  sets,  which  was  particularly 
constraining  for  higher-speed  data  collection  as  well  as  angular  motion. 

Comparing  predicted  and  measured  responses,  the  vehicle  response  close  to  the 
trim  condition  is  accurately  captured,  but  behavior  deviates  farther  away  due  to  the 
complex  dynamics  and  inherent  parameter  dependence  on  vehicle  state.  Dimensionless 
parameters  are  introduced  to  account  for  parameter  deviation  from  the  trim  condition. 
These  dimensionless  parameters  are  learned  with  their  dimensional  counterparts.  The 
dimensionless  parameter  model  formulation  yielded  more  accurate  predictive  capability 
since  the  model  can  correct  for  slight  variations  in  vehicle  velocity  (around  the  trim 
condition).  However,  this  technique  is  not  applicable  to  large  deviations.  To  this  end,  a 
Reynolds  number  analysis  shows  that  the  low-speed  trim  condition  falls  within  the 
laminar-to-transition  flow  regime,  for  which  linear  damping  may  be  a  more  appropriate 
model  (quadratic  damping  is  assumed,  as  motivated  by  vortex  shedding  effects  for 
turbulent  flow).  This  revised  model  is  investigated  for  the  low-speed  trim  condition,  but 
without  significant  improvement  in  model  performance.  As  a  result,  it  is  concluded  that 
operation  in  this  flow  transition  regime  will  likely  require  inclusion  of  both  terms.  The 
original  model  (with  quadratic  damping)  appears  to  be  appropriate  for  the  higher-speed 
operations. 

A  non-parametric  approach  is  also  investigated  based  on  a  computational  NN 
framework.  A  simple  network  is  created  for  both  the  simulator  and  the  SeaBotix 
vLBV300  that  successfully  mapped  the  thruster  force  in  the  surge  direction  to  the  surge 
velocity  u.  Then  a  more  complicated  network  is  created  to  map  the  four  controllable 
DOF  thruster  forces  to  their  respective  body  velocities  and  angular  rate.  This  is 
accomplished  for  both  the  simulator,  but  the  experimental  setup  did  not  allow  sufficient 
data  to  be  captured  for  the  physical  system.  The  potential  of  the  NN  framework  and  its 
ability  to  capture  complex  input-output  relations  inherent  in  modeling  underwater 
vehicles  is  demonstrated.  The  richness  of  the  data  used  to  train  a  NN  is  very  important. 
Much  like  the  concept  of  PE,  a  NN  is  only  as  good  as  the  data  with  which  it  is  trained. 
Then,  the  importance  of  using  a  sufficient  number  of  nodes  and  time  delays  is 
demonstrated,  but  the  pitfall  of  over-fitting  is  also  high-lighted.  Because  this  approach  is 
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non-parametric,  the  resulting  model  does  not  have  the  same  physical  meaning  as  the 
parametric  approaches  presented  so  a  direct  comparison  of  the  identified  parameters  is 
not  possible. 

B.  FUTURE  WORK 

Based  on  the  results  achieved  in  this  work,  as  well  as  the  limitations  identified, 
several  areas  for  future  work  have  been  identified.  First,  the  propulsion  model  is  based 
on  sending  low-level  commands  directly  to  each  thruster.  This  was  beneficial  in  that  it 
allowed  for  control  of  the  vehicle  in  an  additional  DOF,  the  roll  direction.  However,  it  is 
a  limitation  given  the  small  size  of  the  NPS  dive  tank.  Mapping  joystick  commands 
directly  to  the  generated  thrust  will  allow  for  mid-experiment  course  corrections  to  avoid 
collision  or  further  excite  the  vehicle  without  throwing  off  the  force  input  to  the  filters. 

The  assumption  that  the  tether  dynamics  does  not  affect  the  vehicle  dynamics  or 
pose  is  observed  to  be  marginal  at  best  during  the  testing,  particularly  during  operations 
in  the  test  tank.  Further  work  is  required  to  model  the  effects  of  the  tether  to  ensure  the 
more  accurate  control  is  possible. 

CAVR  is  acquiring  an  INS  for  the  SeaBotix  vLBV300  in  the  near  future.  This 
will  overcome  many  of  the  experimental  setup  limitations  encountered  in  this  work.  The 
main  advantage  is  that  the  sensor  will  allow  a  greater  range  of  motions  to  be  tracked,  as 
well  as  ocean-based  operations.  Open  ocean  trials  will  allow  for  much  longer  data  runs 
allowing  the  filters  to  have  more,  and  thanks  to  wave  action,  richer  data  to  use  for 
convergence.  Once  this  has  occurred,  the  approach  and  filters  developed  in  this  work  can 
be  modified  to  use  the  measurements  provided  by  the  INS  in  order  to  improve  the 
accuracy  of  the  filter. 

A  key  assumption  that  allowed  the  use  of  the  classical  RLLS  filter  is  that  the 
parameters  being  estimated  were  not  time  varying.  Because  of  this,  the  RLLS  is 
formulated  to  consider  information  contained  in  the  entire  experiment  (i.e.,  no  old  data  is 
discarded).  Because  the  hydrodynamic  coefficients  were  shown  to  vary  with  vehicle 
speed  and  were  also  shown  to  be  non-scalable  to  that  speed  difference,  a  time-varying 
approach  may  prove  beneficial  and  should  be  investigated.  This  would  not  require 
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substantial  further  work  as  the  classical  RLLS  employed  in  this  work  can  be  readily 
adapted  with  a  “forgetting”  factor  that  allows  it  to  only  consider  data  over  a  certain 
period  of  time.  Other  approaches  for  time-varying  parameters  can  be  investigated. 

Considering  the  Bayesian  filtering  approach:  while  the  system  as  written  is  not 
observable  and  the  extended  Kalman  filter  approach  is  not  usable,  alternative 
mathematical  formulation  of  the  problem  may  produce  the  necessary  condition  for 
observability.  This  will  allow  for  an  investigation  of  not  only  the  extended  Kalman  filter 
approach  but  for  other  statistical  methods  based  on  Bayes’  rule,  such  as  the  unscented 
Kalman  filter  or  a  particle  filter,  which  would  also  prove  useful  when  considering  the 
time-varying  parameter  problem. 
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APPENDIX  A.  6  DOF  MOTION  SIMULATOR  MATLB  CODE 


function  [ydot, etadot, J, out]  =  fen (y, eta, tau) 

%y=  [u;  v;  w;  p;  q;  r  ]  ; 

%tau= [X; Y; Z;K;M;N] ; 

%eta= [x; y; z ; phi ; theta; psi ] ; 

g  =  9.81;  %  Acceleration  due  to  gravity  in  meters/secondA2  (m/sA2) 
m  =  20.9;  %  Mass  in  kilgrams  (kg) 

W  =  m*g;  %  Weight  in  Newtons  (N) 

B  =  W;  %  Measured  Vehicle  Buoyancy  (N) 

%  Moments  of  Inertia  WRT  Origin  at  Half-Length 
I  xx  =  .5587;  %  kg*mA2 
I_yx  =  0; 

I_zx  =  0; 

I_xy  =  0 ; 

I_yy  =  .9531;  %  kg*mA2 
I_yz  =  0; 

I_zy  =  0; 

I  xz=0; 

I  zz  =  .9;  %  kg*mA2 

I=[I  xx  -I  xy  -I  xz;  -I  yx  I  yy  -I  yz;  -I  zx  -I  zy  I  zz]; 

%  Center  of  Buoyancy  WRT  Origin  at  Vehicle  Nose 
x_cb  =  0.00;  %  x-location  (m) 
y_cb  =  0.00;  %  y-location  (m) 
z_cb  =  0.00;  %  z-location  (m) 

%  Center  of  Buoyancy  WRT  Origin  at  Vehicle  Half  Length 
%  half length=L/2 ; 
x_cb  =  0.00; 

y_cb  =  0.00;  %  y-location  (m) 

%  z_cb  =  -5.016E-1;  %  z-location  (m) 
z_cb=- . 05 ; 

%  Center  of  Gravity  WRT  Origin  at  Vehicle  Half  Length 
x  eg  =  0;  %  x-location  (m) 
y_cg  =  0.00;  %  y-location  (m) 
z_cg  =  0;  %  z-location  (m) 

%  Non-Linear  Force  Coefficients 
X  uu  =  -4.56;  %  Cross-flow  Drag  (kg/m) 

X_udot  =  200.4802;  %  Added  Mass  (kg) 

X_u  =  -6.8387; 

X_vdot  =  0; 

X_wdot  =  0; 

X_pdot  =  0; 

X_qdot  =  0; 

X_rdot  =  0; 

X_wq  =0;  %  Added  Mass  Cross-term  (kg/rad) 

X  qq  =  0;  %  Added  Mass  Cross-term  (kg*m/rad) 

X  vr  =0;  %  Added  Mass  Cross-term  (kg/rad) 
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X  rr  =  0;  %  Added  Mass  Cross-term  (kg*m/rad) 

Y_vv  =  -55.006; 

Y  v=0;%  Cross-flow  Drag  (kg/m) 

Y  rr  =  0;  %  Cross-flow  Drag  (kg*m/radA2) 

Y  uv  =  0;  %  Body  Lift  Force  and  Fin  Lift  (kg/m) 

Y_udot  =  0; 

Y  vdot  =  -21.584;  %  Added  Mass  (kg) 

Y_wdot  =  0; 

Y_pdot  =  0; 

Y_qdot  =  0; 

Y_rdot  =0;  %  Added  Mass  (kg*m/rad) 

Y  ur  =0;  %  Added  Mass  Cross-term  and  Fin  Lift  (kg/rad) 

Y_wp  =  0;  %  Added  Mass  Cross-term  (kg/rad) 

Y  pq  =  0;  %  Added  Mass  Cross-term  (kg*m/rad) 

Y  uudr  =  0;  %  Fin  Lift  Force  (kg/ (m*rad) ) 

Z  ww  =  -67.8358; 

Z  w=0;%  Cross-flow  Drag  (kg/m) 

Z  qq  =  0;  %  Cross-flow  Drag  (kg*m/rad) 

Z  uw  =  0;  %  Body  Lift  Force  and  Fin  Lift  (kg/m) 

Z_udot  =  0; 

Z_vdot  =  0; 

Z  wdot  =  -21.3775;  %  Added  Mass  (kg) 

Z_pdot  =  0; 

Z_qdot  =0;  %  Added  Mass  (kg*m/rad) 

Z_rdot  =  0; 

Z  uq  =  0;  %  Added  Mass  Cross-term  and  Fin  Lift  (kg/rad) 

Z  vp  =0;  %  Added  Mass  Cross-term  (kg/rad) 

Z_rp  =0;  %  Added  Mass  Cross-term  (kg/rad) 

Z  uuds  =  0;  %  Fin  Lift  Force  (kg/ (m*rad) ) 

K_pp  =  -103.335; 

K  p=0;%  Rolling  Resistance  ( kg*mA2 /radA2 ) 

K_udot  =  0; 

K_vdot  =  0; 

K_wdot  =  0; 

K  pdot  =  -23.890;  %  Added  Mass  (kg*mA2/rad) 

K_qdot  =  0; 

K_rdot  =  0; 

M  ww  =  0;  %  Cross-flow  Drag  (kg) 

M~qq  =  -129.99; 

M  q=0;%  Cross-flow  Drag  ( kg*mA2 /radA2 ) 

M  uw  =0;  %  Body  and  Fin  Lift  and  Munk  Moment  (kg) 

M_udot  =  0; 

M_vdot  =  0; 

M_wdot  =0;  %  Added  Mass  (kg*m) 

M_pdot  =  0; 

M  qdot  =  -46.726;  %  Added  Mass  (kg*mA2/rad) 

M_rdot  =  0; 

M  uq  =  0;  %  Added  Mass  Cross-term  and  Fin  Lift  (kg*m/rad) 
M_vp  =0;  %  Added  Mass  Cross-term  (kg*m/rad) 

M  rp  =  0;  %  Added  Mass  Cross-term  ( kg*mA2 /radA2 ) 

M  uuds  =0;  %  Fin  Lift  Moment  (kg/rad) 

N  vv  =  0;  %  Cross-flow  Drag  (kg) 

N~rr  =  -.8814; 

N  r=0;%  Cross-flow  Drag  ( kg*mA2 /radA2 ) 
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N  uv  =  0;  %  Body  and  Fin  Lift  and  Munk  Moment  (kg) 

N_udot  =  0; 

N  vdot  =0;  %  Added  Mass  (kg*m) 

N  wdot  =  0; 

N_pdot  =  0; 

N_qdot  =  0; 

N  rdot  =  -.4202;  %  Added  Mass  (kg*mA2/rad) 

N  ur  =  0;  %  Added  Mass  Cross-term  and  Fin  Lift  (kg*m/rad) 

N  wp  =  0;  %  Added  Mass  Cross-term  (kg*m/rad) 

N  pq  =  0;  %  Added  Mass  Cross-term  ( kg*mA2 /radA2 ) 

N  uudr  =0;  %  Fin  Lift  Moment  (kg/rad) 

tauf=tau; 

%%  Mass  Matrices 

Ma  =  -1*[X  udot  X  vdot  X  wdot  X  pdot  X  qdot  X  rdot; . . . 
Y_udot  Y_vdot  Y_wdot  Y_pdot  Y_qdot  Y_rdot; . . . 

Z  udot  Z  vdot  Z  wdot  Z  pdot  Z  qdot  Z  rdot; . . . 

K  udot  K  vdot  K  wdot  K  pdot  K  qdot  K  rdot; . . . 

M  udot  M  vdot  M  wdot  M  pdot  M  qdot  M  rdot; . . . 

N  udot  N  vdot  N  wdot  N  pdot  N  qdot  N  rdot] ; 

Mrb  =  [m  0  0  0  m*z  eg  -m*y_cg; . . . 

0  m  0  -m*z  eg  0  m*x  eg;  .  . . 

0  0m  m*y  eg  -m*x  eg  0;  .  .  . 

0  -m*z  eg  m*y  eg  I  xx  -I  xy  -I  xz; . . . 
m*z  eg  0  -m*x  eg  -I  yx  I  yy  -I  yz;  .  .  . 

-m*y_cg  m*x  eg  0  -I  zx  -I  zy  I  zz]  ; 

M  =  Ma+Mrb; 

%%  Coriolis  and  Centripetal  Matrix  Calculation 
%y= [u; v; w; p; q; r ]  ; 

%tau= [X; Y; Z;K;M;N] ; 

%eta= [x; y; z ; phi ; theta; psi ] ; 
u=y ( 1 ) ; 
v=y (2) ; 
w=y ( 3 ) ; 

P=y(4) ; 
q=y ( 5 ) ; 
r=y (6) ; 

al  =  X  udot*u+X  vdot*v+X  wdot*w+X  pdot*p+X  qdot*q+X  rdot*r 

a2  =  X  vdot*u+Y  vdot*v+Y  wdot*w+Y  pdot*p+Y  qdot*q+Y  rdot*r 

a3  =  X  wdot*u+Y  wdot*v+Z  wdot*w+Z  pdot*p+Z  qdot*q+Z  rdot*r 

bl  =  X  pdot*u+Y  pdot*v+Z  pdot*w+K  pdot*p+K  qdot*q+K  rdot*r 

b2  =  X  qdot*u+Y  qdot*v+Z  qdot*w+K  qdot*p+M  qdot*q+M  rdot*r 

b3  =  X  rdot*u+Y  rdot*v+Z  rdot*w+K  rdot*p+M  rdot*q+N  rdot*r 


[0 

0 

0 

0 

-a3 

a2  ; 

0 

0 

0 

a3 

0 

-al; 

0 

0 

0 

-a2 

al 

0;  . 

0 

-a3 

a2 

0 

-b3 

b2  ; 

a3 

0 

-al 

b3 

0 

-bl; 

-a2 

al 

0 

-b2 

bl 

0]  ; 

%Fossen's  Simplified  for  UUV 
rbl  =  zeros ( 3 ) ; 

rb2  =  [-m* (y_cg*q+z  cg*r)  m* (y_cg*p+w)  m* (z  cg*p-v) ; . . . 
m* (x  cg*q-w)  -m* (z  cg*r+x  cg*p)  m* (z  cg*q+u) ; . . . 
m*  (x  cg*r+v)  m* (y  cg*r-u)  -m* (x  cg*p+y  cg*q) ] ; 
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rb3  =  [m* (y  cg*q+z  cg*r)  -m* (x  cg*q-w)  -m* (x  cg*r+v) ; . . . 
-m* (y  cg*p+w)  m* (z  cg*r+x  cg*p)  -m* (y  cg*r-u) ; . . . 
-m*(z  cg*p-v)  -m*(z  cg*q+u)  m* (x  cg*p+y  cg*q) ] ; 
rb4  =  [0  -I_yz*q-I  xz*p+I  zz*r  I  yz*r+I  xy*p-I  yy*q; . . . 

I  yz*q+I  xz*p-I  zz*r  0  -I  xz*r-I  xy*q+I  xx*p; . . . 

-I  yz*r-I  xy*p+I  yy*q  I  xz*r+I  xy*q-I  xx*p  0]; 

Crb  =  [rbl  rb3;  rb2  rb4]; 


C  =  Ca+Crb; 

%%  Damping  Matrix 
dl  =  [X_u  00000;. 
0  Y_v  0  0  0  0 ; . . . 
0  0  Z_w  0  0  0 ;  . . . 
0  0  0  K_p  0  0 ;  .  .  . 
0000  M_q  0 ; . . . 
00000  N_r] ; 
d2  =  [X_uu*abs (u)  0  0 
0  Y  vv*abs (v)  0  0 
0  0  Z  ww*abs (w)  0 
000  K_pp*abs (p) 
0  0  M  ww*abs (w)  0 
0  N  vv*abs (v)  0  0 
D  =-dl ; 


0  0  0 ;  .  .  . 

0  Y_rr*abs (r) ; . . . 
Z_qq*abs (q)  0; . . . 
0  0;  .  .  . 

M_qq*abs (q)  0; . . . 
0  N  rr*abs (r) ] ; 


out— [D ( 1 , 1 ) ,  dl ( 1 , 1 ) , d2 ( 1 , 1 ) , u] ; 


%%  Restoring  Force  Matrix 
%eta= [x; y; z;phi; theta;psi] ; 
phi=eta ( 4 ) ; 
theta=eta ( 5 ) ; 
psi=eta ( 6) ; 

gn  =  [ (W-B) *sin (theta) ;  .  .  . 

- (W-B) *cos (theta) *sin (phi) ; .  .  . 

- (W-B) *cos (theta) *cos (phi) ; . . . 

-(y  cg*W-y  cb*B) *cos (theta) *cos (phi) + (z  cg*W- 
z_cb*B) *cos (theta) *sin (phi) ; . . . 

(z  cg*W-z  cb*B) *sin (theta) + (x  cg*W-x  cb*B) *cos (theta) *cos (phi) . 
-(x  cg*W-x  cb*B) *cos (theta) *sin (phi) - (y  cg*W-y  cb*B) *sin  (theta) ] ; 
gnd=gn; 

%%  Y  dot  calculation 
ydot  =  M\(tauf-C*y-D*y-gn); 

Jl= [cos (psi) *cos (theta)  -sin (psi) *cos (phi) +cos (psi) *sin (theta) *sin (phi) 
sin (psi) *sin (phi) +cos (psi) *cos (phi) *sin (theta) ; . . . 

sin (psi) *cos (theta)  cos  (psi) *cos (phi) +sin (phi) *sin (theta) * sin (psi) 
-cos (psi) *sin (phi) +sin (theta) *sin (psi) *cos (phi) ; . . . 

-sin (theta)  cos (theta) *sin (phi)  cos (theta) *cos (phi) ] ; 

J2=  [1  sin (phi) *tan (theta)  cos (phi ) *tan (theta) . 

0  cos (phi)  -sin  (phi) ; . . . 

0  sin (phi) /cos (theta)  cos (phi) /cos (theta) ] ; 

J=[J1,  zeros (3);  zeros ( 3 ) , J2 ] ; 


etadot  =  J*y; 
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APPENDIX  B.  SIMULINK  DIAGRAMS 


Figure  28.  RLLS  estimator 


0— *9i 


surge  velocity 


it 


Band-Limited  White  Noise 


-*E>- 


£>• 


lit 


Regressor 


□ 


Matrix 

Multiply 


To  Workspace 


n 


/* 


Scope 


Figure  29.  Gradient  estimator 
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Figure  30.  Simulator 
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